Introduction

Here we will be running the analysis for the hemp core microbiome project. Anka hemp plants were grown in 6 fields around Ithaca and Geneva NY. We collected bulk soil, rhizosphere soil, root tissue, leaf washes, and flower washes from 5 plants per field. We extracted DNA from these samples and amplified and sequenced the V4 region of the 16S rRNA gene (bacterial community) as well as the ITS1 region (fungal community). The goal of this study is to identify members of the core microbiome of hemp.

Calling packages

# Packages needed for data handeling
library(dplyr)
library(tibble)
library(tidyr)
library(phyloseq)
library(kableExtra)
library(multcomp)

# Packages needed for analysis
library(vegan)
library(DESeq2)
library(ape)
library(knitr)
library(agricolae)

# Packages needed for plotting
library(ggplot2)

# Setting seed similar to R versions < 3.6
RNGkind(sample.kind = "Rounding")

Importing data

All data is combined in phyloseq objects. First we need to do a bit of filtering

# Import field data including crop yield
field.data = read.table("/home/sam/notebooks/hemp_microbiome/data/hemp_field_data.txt", stringsAsFactors = F, sep="\t", header=T)

# Import bacterial/archeal phyloseq object created in the sequence analysis pipeline
physeq = readRDS("/home/sam/notebooks/hemp_microbiome/data/16S_OTUs/physeq.hemp16S.rds")

# Get the metadata table for later use
metadata = data.frame(sample_data(physeq), stringsAsFactors = F)

# Extract just OTUs classified as bacteria
physeq16S.bact = subset_taxa(physeq, Rank1=="D_0__Bacteria")
physeq = NULL # Removing the original phyloseq to save memory

# I'll also import the phylogenetic tree for the bacterial OTUs and add it to the bacterial phyloseq object
bact.tree = read.tree("/home/sam/notebooks/hemp_microbiome/data/16S_OTUs/hemp.bact.tre")
phy_tree(physeq16S.bact) = bact.tree

# Import the fungal phyloseq object
physeqITS.fungi = readRDS("/home/sam/notebooks/hemp_microbiome/data/ITS_OTUs/physeq.hempITS.rds")

# Remove any OTUs not classified as Fungi
physeqITS.fungi = subset_taxa(physeqITS.fungi, Kingdom == "Fungi")


# Remove unsequenced samples. These were samples that were not sequenced due to low PCR product. They were originally left in the metadata file for simplicity sake.
physeq16S.bact = subset_samples(physeq16S.bact, Sequenced_16S == "TRUE")
physeq16S.bact = prune_taxa(taxa_sums(physeq16S.bact) > 0, physeq16S.bact)
physeqITS.fungi = subset_samples(physeqITS.fungi, Sequenced_ITS == "TRUE")
physeqITS.fungi = prune_taxa(taxa_sums(physeqITS.fungi) > 0, physeqITS.fungi)

# Remove sample F5-1 because it does not look right. In the ordination using flowers, it falls way far away. This is probably due to contamination.
physeq16S.bact = subset_samples(physeq16S.bact, sample_ID != "F5-1")

Now remove any samples that are poorly sequenced (<2000 reads for bacteria and <1000 reads for fungi).

## Find samples that have OTU totals greater than a threshold (2000 for bacteria and 1000 for fungi) and subset only these
physeq16S.bact = prune_samples(sample_sums(physeq16S.bact) >= 2000, physeq16S.bact)
physeq16S.bact = prune_taxa(taxa_sums(physeq16S.bact) > 0, physeq16S.bact)

physeqITS.fungi = prune_samples(sample_sums(physeqITS.fungi) >= 1000, physeqITS.fungi)
physeqITS.fungi = prune_taxa(taxa_sums(physeqITS.fungi) > 0, physeqITS.fungi)

Also need to remove some extra sample types that were not used in this analysis and root tissues from the fungal dataset because they were primarily hemp ITS sequences.

## Remove the Male/Female flower samples. Also remove the root tissue samples from the fungal dataset since these were overrun with hemp ITS sequences.
physeq16S.bact = subset_samples(physeq16S.bact, !(Sample_type == "male_flowers" | Sample_type == "female_flowers"))
physeq16S.bact = prune_taxa(taxa_sums(physeq16S.bact) > 0, physeq16S.bact)
physeqITS.fungi = subset_samples(physeqITS.fungi, !(Sample_type == "male_flowers" | Sample_type == "female_flowers" | Sample_type == "root_tissue"))
physeqITS.fungi = prune_taxa(taxa_sums(physeqITS.fungi) > 0, physeqITS.fungi)

## Also remove the water controls
physeq16S.bact = subset_samples(physeq16S.bact, !(Sample_type == "water_control"))
physeq16S.bact = prune_taxa(taxa_sums(physeq16S.bact) > 0, physeq16S.bact)
physeqITS.fungi = subset_samples(physeqITS.fungi, !(Sample_type == "water_control"))
physeqITS.fungi = prune_taxa(taxa_sums(physeqITS.fungi) > 0, physeqITS.fungi)

The read depth of our samples varies quite a bit. To normalize this lets use rarefaction.

## Rarify samples. This samples the OTUs for all samples to an even depth to account for large differences in sample sizes. Sampling is based 
set.seed(72)
physeq16S.bact.rare = rarefy_even_depth(physeq16S.bact)
physeq16S.bact.rare = prune_taxa(taxa_sums(physeq16S.bact.rare) > 0, physeq16S.bact.rare)
set.seed(72)
physeqITS.fungi.rare = rarefy_even_depth(physeqITS.fungi)
physeqITS.fungi.rare = prune_taxa(taxa_sums(physeqITS.fungi.rare) > 0, physeqITS.fungi.rare)

## Normalize the OTU count data to the total for that sample. This gives you relative abundance for each OTU.
physeq16S.bact.rare = transform_sample_counts(physeq16S.bact.rare, function(x) x/sum(x))
physeqITS.fungi.rare = transform_sample_counts(physeqITS.fungi.rare, function(x) x/sum(x))

Seting colors, shapes, and figure functions

I want to have a consistent color and shape scheme through all figures so I will set that here. Colors come from https://personal.sron.nl/~pault/data/colourschemes.pdf and should be colorblind friendly.

# Set colors
## Overlapping
type.col = c("Bulk soil" = "#1965B0", "Flowers" = "#7BAFDE", "Leaves" = "#4EB265", "Rhizosphere soil" = "#F7F056", "Root tissue" = "#DC050C")
field.col = c("Crittenden North" = "#1965B0", "McGowan Late" = "#7BAFDE", "East Ithaca" = "#4EB265", "McGowan" = "#CAE0AB", "Research North" = "#F7F056", "Freeville" = "#DC050C")

## Non-overlapping. We will not use these
#type.col = c("Bulk soil" = "#882E72", "Flowers" = "#1965B0", "Leaves" = "#5289C7", "Rhizosphere soil" = "#7BAFDE", "Root tissue" = "#4EB265")
#field.col = c("Crittenden North" = "#CAE0AB", "McGowan Late" = "#F7F056", "East Ithaca" = "#F4A736", "McGowan" = "#E8601C", "Research North" = "#DC050C", "Freeville" = "#72190E")

## What does this look like
pie(rep(1,5), labels=names(type.col), col=type.col)

pie(rep(1,6), labels=names(field.col), col=field.col)

# Set Shapes
field.shape = c("Crittenden North" = 21, "McGowan Late" = 22, "East Ithaca" = 23, "McGowan" = 24, "Research North" = 25, "Freeville" = 8)

# This is just a function for getting the legends for graphs. This will be used when making plots for publication.
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

# These are functions for reordering the x axis within faceted figures
reorder_within <- function(x, by, within, fun = mean, sep = "___", ...) {
  new_x <- paste(x, within, sep = sep)
  stats::reorder(new_x, by, FUN = fun)
}
scale_x_reordered <- function(..., sep = "___") {
  reg <- paste0(sep, ".+$")
  ggplot2::scale_x_discrete(labels = function(x) gsub(reg, "", x), ...)
}

Part 1: Community composition and ordination analysis

In this analysis I will compare the overall bacterial and fungal community compositions between sample types (i.e. plant parts) and fields.

Inital stats

First lets see how many OTUs there are and how many phyla are found

print(paste("There are", ntaxa(physeq16S.bact.rare), "bacterial OTUS"))
## [1] "There are 8913 bacterial OTUS"
print(paste("There are", ntaxa(physeqITS.fungi.rare), "fungal OTUS"))
## [1] "There are 982 fungal OTUS"
print(paste("There are", length(na.omit(c(unique(tax_table(physeq16S.bact.rare)[,2])))), "bacterial Phyla"))
## [1] "There are 39 bacterial Phyla"
print(paste("There are", length(na.omit(c(unique(tax_table(physeqITS.fungi.rare)[,2])))), "fungal Phyla"))
## [1] "There are 9 fungal Phyla"

Alpha diversity

Does OTU evenness vary across plant compartments?

new_sample_types = data.frame(Sample_type = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"),
                              Sample_type_new = c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))

bact.alphadiv = full_join(tibble::rownames_to_column(data.frame(Shannon = diversity(t(otu_table(physeq16S.bact.rare)), index = "shannon")), var="sample_ID"),
                         tibble::rownames_to_column(data.frame(Richness = specnumber(t(otu_table(physeq16S.bact.rare)))), var="sample_ID"), by = "sample_ID") %>%
  mutate(Evenness = Shannon/log(Richness)) %>%
  left_join(data.frame(sample_data(physeq16S.bact.rare)), by = "sample_ID") %>%
  left_join(new_sample_types, by="Sample_type") %>%
  mutate(Sample_type_new = factor(Sample_type_new, 
                                  levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers")),
         Sample_type = factor(Sample_type, 
                              levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers")),
         Field = factor(Field), organism = "Bacteria") %>%
  dplyr::select(organism, sample_ID, Shannon, Richness, Evenness, Field, Sample_type, Sample_type_new) %>%
  tidyr::gather(key = "Index", value = "diversity", -organism, -sample_ID, -Field, -Sample_type, -Sample_type_new)

fungi.alphadiv = full_join(tibble::rownames_to_column(data.frame(Shannon = diversity(t(otu_table(physeqITS.fungi.rare)), index = "shannon")), var="sample_ID"),
                         tibble::rownames_to_column(data.frame(Richness = specnumber(t(otu_table(physeqITS.fungi.rare)))), var="sample_ID"), by = "sample_ID") %>%
  mutate(Evenness = Shannon/log(Richness)) %>%
  left_join(data.frame(sample_data(physeqITS.fungi.rare)), by = "sample_ID") %>%
  left_join(new_sample_types, by="Sample_type") %>%
  mutate(Sample_type_new = factor(Sample_type_new, 
                                  levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers")),
         Sample_type = factor(Sample_type, 
                              levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers")),
         Field = factor(Field), organism = "Fungi") %>%
  dplyr::select(organism, sample_ID, Shannon, Richness, Evenness, Field, Sample_type, Sample_type_new) %>%
  tidyr::gather(key = "Index", value = "diversity", -organism, -sample_ID, -Field, -Sample_type, -Sample_type_new)
  
alphadiv = rbind(bact.alphadiv, fungi.alphadiv)

div.aov.list = list()
div.cld.df = data.frame()

for (org in c("Bacteria", "Fungi")){
  div_index = "Evenness"
  div.model = lm(diversity ~ Sample_type, data=filter(alphadiv, Index == div_index, organism == org))
  div.aov = anova(div.model)
  div.aov.list[[paste(org, div_index, sep="_")]] = div.aov
  div.posthoc = glht(div.model, linfct = mcp(Sample_type = "Tukey"))
  div.cld = data.frame(Index = div_index, sig_group = cld(div.posthoc)$mcletters$Letters) %>%
    tibble::rownames_to_column(var="Sample_type") %>%
    left_join(new_sample_types, by="Sample_type") %>%
    mutate(Sample_type_new = factor(Sample_type_new, 
                                    levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers")),
           organism = org)
  div.cld.df = rbind(div.cld.df, div.cld)
}

div.aov.list
## $Bacteria_Evenness
## Analysis of Variance Table
## 
## Response: diversity
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Sample_type   4 2.7370 0.68425  60.657 < 2.2e-16 ***
## Residuals   130 1.4665 0.01128                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Fungi_Evenness
## Analysis of Variance Table
## 
## Response: diversity
##             Df  Sum Sq  Mean Sq F value    Pr(>F)    
## Sample_type  3 0.85717 0.285722  53.059 < 2.2e-16 ***
## Residuals   88 0.47388 0.005385                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alphadiv = rbind(alphadiv, 
                 data.frame(organism = "Fungi", sample_ID = NA, Field = NA, Sample_type = NA,
                            Sample_type_new = "Root tissue", Index = "Evenness", diversity=NA))

bact_alpha.plot = ggplot(data=filter(alphadiv, Index == "Evenness", organism == "Bacteria"), aes(x=Sample_type_new, y=diversity)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.3) +
  geom_text(data=filter(div.cld.df, Index == "Evenness", organism == "Bacteria"), aes(x=Sample_type_new, label=sig_group), y=1, size=4) +
  lims(y=c(0, 1)) +
  labs(x="Compartment", y="Pielou's Evenness") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        text = element_text(size=16))

fungi_alpha.plot = ggplot(data=filter(alphadiv, Index == "Evenness", organism == "Fungi"), aes(x=Sample_type_new, y=diversity)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha=0.3) +
  geom_text(data=filter(div.cld.df, Index == "Evenness", organism == "Fungi"), aes(x=Sample_type_new, label=sig_group), y=1, size=4) +
  lims(y=c(0, 1)) +
  labs(x="Compartment", y="Pielou's Evenness") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        text = element_text(size=16))


# Use cowplot::plot_grid to plot both together. Remove the legend from the plots and add in the full legend
alpha.plot = cowplot::plot_grid(bact_alpha.plot + theme(axis.title = element_blank(), axis.text.x = element_blank(),
                                                        plot.margin = margin(7, 0, 0, 5, "mm")), 
                                fungi_alpha.plot + theme(axis.title = element_blank(),
                                                         plot.margin = margin(7, 0, 0, 5, "mm")), 
                                align = "v", rel_heights = c(0.7, 1), ncol=1, labels = c(" a", " b"))

y.grob <- grid::textGrob("                     Pielou's Evenness", gp=grid::gpar(fontsize=14), rot=90)
x.grob <- grid::textGrob("          Compartment", gp=grid::gpar(fontsize=14))

alpha.plot = gridExtra::arrangeGrob(alpha.plot, left = y.grob, bottom = x.grob)

gridExtra::grid.arrange(alpha.plot)

## Saving figure for publication
options(bitmapType='cairo')
#ggsave(plot = alpha.plot, filename = "/home/sam/notebooks/hemp_microbiome/figures/evenness_plot.tiff",
#       device = "tiff", width=3, height=6, units="in")

Taxanomic makeup

Now lets visuallize the taxonomic makeup of the communities at the phylum or class level. Simply, I will calculate the proportion of reads from across all plants/samples that are assigned to each phylum or class. I’ll do this separately for each sample type so that we can compare the rough composition between types.

Bacteria

bact.phyla.sum = data.frame()

for (samtype in c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers")){
  sub.physeq = subset_samples(physeq16S.bact.rare, Sample_type == as.character(samtype))
  sub.physeq = prune_taxa(taxa_sums(sub.physeq) > 0, sub.physeq)
  sub.otu.table = otu_table(sub.physeq)
  sub.tax.table = data.frame(tax_table(sub.physeq), stringsAsFactors = FALSE) %>%
    rownames_to_column(var="OTU")
  sub.otu.df = data.frame(sub.otu.table) %>%
    rownames_to_column(var="OTU") %>%
    gather(key="sample", value="count", -OTU) %>% 
    filter(count > 0) %>%
    inner_join(sub.tax.table, by="OTU")
  
  sub.phyla.sum = sub.otu.df %>%
    mutate(taxa = gsub("D_.__", "", ifelse(Rank2 == "D_1__Proteobacteria", as.character(Rank3), as.character(Rank2)))) %>%
    group_by(taxa, sample) %>%
    summarize(RA = sum(count)*100) %>%
    as.data.frame() %>%
    tidyr::spread(key=sample, value=RA) %>%
    tidyr::gather(key=sample, value=RA, -taxa) %>%
    mutate(RA = ifelse(is.na(RA), 0, RA)) %>%
    mutate(Sample_type = samtype)
  
  bact.phyla.sum = rbind(bact.phyla.sum, sub.phyla.sum)
}

bact.phyla.sum.mean = bact.phyla.sum %>%
  group_by(taxa, Sample_type) %>%
  summarize(mean_RA = mean(RA)) %>%
  as.data.frame
  
bact.high.taxa = unique(bact.phyla.sum.mean[bact.phyla.sum.mean$mean_RA >= 1,]$taxa)
bact.taxa.colors = c("#D1BBD7", "#AE76A3", "#882E72", "#1965B0", "#5289C7", "#7BAFDE", "#4EB265", "#90C987", "#CAE0AB", "#F7F056", "#F4A736", "#E8601C", "#DC050C", "grey")
names(bact.taxa.colors) = c(bact.high.taxa, "less than 1% each")

new_sample_types = data.frame(Sample_type = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"),
                              Sample_type_new = c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))

bact.phyla.sum.high = bact.phyla.sum %>%
  mutate(taxa = ifelse(taxa %in% bact.high.taxa, taxa, "less than 1% each")) %>%
  group_by(taxa, sample, Sample_type) %>%
  summarize(RA = sum(RA)) %>%
  as.data.frame %>%
  group_by(taxa, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA, na.rm = TRUE),
            n_samples = n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  inner_join(new_sample_types)
bact.phyla.sum.high$taxa = factor(bact.phyla.sum.high$taxa, levels=c(sort(bact.high.taxa), "less than 1% each"))
bact.phyla.sum.high$Sample_type_new = factor(bact.phyla.sum.high$Sample_type_new, levels=c("Flowers", "Leaves", "Root tissue", "Rhizosphere soil", "Bulk soil"))


bact.comm.plot = ggplot(data=bact.phyla.sum.high, aes(x=taxa, y=mean_RA, fill=taxa)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=mean_RA-SE_RA, ymax=mean_RA+SE_RA), width = 0.5, size=0.25) +
  scale_fill_manual(values=bact.taxa.colors) +
  labs(x="Phylum/Class", y="Relative abundance (% of reads)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  facet_wrap(~Sample_type_new, nrow = 5)
bact.comm.plot

Here is the same data but in table format.

bact.phyla.sum.high.df = bact.phyla.sum.high %>%
  dplyr::select(taxa, Sample_type, mean_RA) %>%
  spread(key=Sample_type, value=mean_RA)

kable(bact.phyla.sum.high.df, "html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
taxa bulk_soil flowers leaves rhizosphere_soil root_tissue
Acidobacteria 12.623315 3.6161776 1.0073638 13.650957 3.8930890
Actinobacteria 26.910389 9.5463918 15.7069219 27.020619 29.7991600
Alphaproteobacteria 10.517050 10.0570975 44.2695140 10.044183 13.0798015
Bacteroidetes 2.812054 2.2727994 17.7407953 2.543446 6.6498664
Betaproteobacteria 7.858842 5.8445678 4.6318115 8.839470 17.9488354
Chloroflexi 8.867565 2.0269627 0.8424153 8.243004 3.1660939
Deltaproteobacteria 4.301348 1.5384615 1.3357879 4.092784 2.0771287
Firmicutes 2.126883 22.6724822 6.0868925 1.860088 4.6170294
Gammaproteobacteria 5.162569 38.1601903 6.0250368 5.019146 14.1825124
Gemmatimonadetes 5.960349 1.1102300 0.4226804 5.329897 0.8705613
Nitrospirae 1.595559 0.2442506 0.0810015 1.863034 0.2794960
Planctomycetes 5.053132 1.3291039 0.4064801 5.656848 1.2523864
Verrucomicrobia 3.395718 0.7359239 0.3490427 3.836524 1.1347843
less than 1% each 2.815226 0.8453608 1.0942563 2.000000 1.0492554

Using nonparametric Kruskal test with Dunn posthoc test, lets see if the relative abundances of these phyla differ across plant compartments.

high.taxa = unique(bact.phyla.sum.mean[bact.phyla.sum.mean$mean_RA >= 1,]$taxa)
bact.posthoc.res = data.frame()
for (subtaxa in high.taxa){
  taxa.model = kruskal.test(RA~Sample_type, data=filter(bact.phyla.sum, taxa == subtaxa))
  if (p.adjust(taxa.model$p.value, n=length(high.taxa), method="bonferroni") < 0.05){
    posthoc.res = FSA::dunnTest(RA~Sample_type, data=filter(bact.phyla.sum, taxa == subtaxa))
    posthoc.res = posthoc.res$res
    posthoc.cld = rcompanion::cldList(comparison = posthoc.res$Comparison,
                                      p.value    = posthoc.res$P.adj,
                                      threshold  = 0.05)
  } else{
    posthoc.cld = data.frame(Group = NA, Letter = NA)
  }
  sub.bact.posthoc.res = data.frame(taxa = subtaxa, X2 = taxa.model$statistic,
                                    Pvalue = taxa.model$p.value,
                                    padj = p.adjust(taxa.model$p.value, n=length(high.taxa), method="bonferroni"),
                                    Sample_type = posthoc.cld$Group,
                                    group = posthoc.cld$Letter,
                                    stringsAsFactors = FALSE)
  bact.posthoc.res = rbind(bact.posthoc.res, sub.bact.posthoc.res)
}
bact.posthoc.res$Sample_type = factor(bact.posthoc.res$Sample_type, levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"))
#write.table(bact.posthoc.res, "/home/sam/notebooks/hemp_microbiome/data/bact_phyla_across_compartments.txt", quote=FALSE, row.names = FALSE)

kable(bact.posthoc.res %>%
        arrange(taxa, Sample_type), "html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
taxa X2 Pvalue padj Sample_type group
Acidobacteria 92.86427 0.0e+00 0.00e+00 bulk_soil a
Acidobacteria 92.86427 0.0e+00 0.00e+00 rhizosphere_soil a
Acidobacteria 92.86427 0.0e+00 0.00e+00 root_tissue c
Acidobacteria 92.86427 0.0e+00 0.00e+00 leaves b
Acidobacteria 92.86427 0.0e+00 0.00e+00 flowers bc
Actinobacteria 68.39317 0.0e+00 0.00e+00 bulk_soil a
Actinobacteria 68.39317 0.0e+00 0.00e+00 rhizosphere_soil a
Actinobacteria 68.39317 0.0e+00 0.00e+00 root_tissue a
Actinobacteria 68.39317 0.0e+00 0.00e+00 leaves b
Actinobacteria 68.39317 0.0e+00 0.00e+00 flowers b
Alphaproteobacteria 68.52089 0.0e+00 0.00e+00 bulk_soil a
Alphaproteobacteria 68.52089 0.0e+00 0.00e+00 rhizosphere_soil a
Alphaproteobacteria 68.52089 0.0e+00 0.00e+00 root_tissue a
Alphaproteobacteria 68.52089 0.0e+00 0.00e+00 leaves b
Alphaproteobacteria 68.52089 0.0e+00 0.00e+00 flowers a
Bacteroidetes 81.51694 0.0e+00 0.00e+00 bulk_soil a
Bacteroidetes 81.51694 0.0e+00 0.00e+00 rhizosphere_soil a
Bacteroidetes 81.51694 0.0e+00 0.00e+00 root_tissue b
Bacteroidetes 81.51694 0.0e+00 0.00e+00 leaves b
Bacteroidetes 81.51694 0.0e+00 0.00e+00 flowers a
Betaproteobacteria 69.61585 0.0e+00 0.00e+00 bulk_soil ab
Betaproteobacteria 69.61585 0.0e+00 0.00e+00 rhizosphere_soil b
Betaproteobacteria 69.61585 0.0e+00 0.00e+00 root_tissue d
Betaproteobacteria 69.61585 0.0e+00 0.00e+00 leaves c
Betaproteobacteria 69.61585 0.0e+00 0.00e+00 flowers ac
Chloroflexi 96.14107 0.0e+00 0.00e+00 bulk_soil a
Chloroflexi 96.14107 0.0e+00 0.00e+00 rhizosphere_soil a
Chloroflexi 96.14107 0.0e+00 0.00e+00 root_tissue c
Chloroflexi 96.14107 0.0e+00 0.00e+00 leaves b
Chloroflexi 96.14107 0.0e+00 0.00e+00 flowers bc
Deltaproteobacteria 61.18579 0.0e+00 0.00e+00 bulk_soil a
Deltaproteobacteria 61.18579 0.0e+00 0.00e+00 rhizosphere_soil a
Deltaproteobacteria 61.18579 0.0e+00 0.00e+00 root_tissue b
Deltaproteobacteria 61.18579 0.0e+00 0.00e+00 leaves b
Deltaproteobacteria 61.18579 0.0e+00 0.00e+00 flowers b
Firmicutes 32.17206 1.8e-06 2.29e-05 bulk_soil a
Firmicutes 32.17206 1.8e-06 2.29e-05 rhizosphere_soil a
Firmicutes 32.17206 1.8e-06 2.29e-05 root_tissue a
Firmicutes 32.17206 1.8e-06 2.29e-05 leaves a
Firmicutes 32.17206 1.8e-06 2.29e-05 flowers b
Gammaproteobacteria 61.02114 0.0e+00 0.00e+00 bulk_soil a
Gammaproteobacteria 61.02114 0.0e+00 0.00e+00 rhizosphere_soil a
Gammaproteobacteria 61.02114 0.0e+00 0.00e+00 root_tissue b
Gammaproteobacteria 61.02114 0.0e+00 0.00e+00 leaves a
Gammaproteobacteria 61.02114 0.0e+00 0.00e+00 flowers b
Gemmatimonadetes 95.69993 0.0e+00 0.00e+00 bulk_soil a
Gemmatimonadetes 95.69993 0.0e+00 0.00e+00 rhizosphere_soil a
Gemmatimonadetes 95.69993 0.0e+00 0.00e+00 root_tissue b
Gemmatimonadetes 95.69993 0.0e+00 0.00e+00 leaves b
Gemmatimonadetes 95.69993 0.0e+00 0.00e+00 flowers b
Nitrospirae 89.15074 0.0e+00 0.00e+00 bulk_soil a
Nitrospirae 89.15074 0.0e+00 0.00e+00 rhizosphere_soil a
Nitrospirae 89.15074 0.0e+00 0.00e+00 root_tissue b
Nitrospirae 89.15074 0.0e+00 0.00e+00 leaves b
Nitrospirae 89.15074 0.0e+00 0.00e+00 flowers b
Planctomycetes 91.24467 0.0e+00 0.00e+00 bulk_soil a
Planctomycetes 91.24467 0.0e+00 0.00e+00 rhizosphere_soil a
Planctomycetes 91.24467 0.0e+00 0.00e+00 root_tissue c
Planctomycetes 91.24467 0.0e+00 0.00e+00 leaves b
Planctomycetes 91.24467 0.0e+00 0.00e+00 flowers bc
Verrucomicrobia 91.51265 0.0e+00 0.00e+00 bulk_soil a
Verrucomicrobia 91.51265 0.0e+00 0.00e+00 rhizosphere_soil a
Verrucomicrobia 91.51265 0.0e+00 0.00e+00 root_tissue c
Verrucomicrobia 91.51265 0.0e+00 0.00e+00 leaves b
Verrucomicrobia 91.51265 0.0e+00 0.00e+00 flowers bc

Plot the results.

bact.phyla.high = bact.phyla.sum %>%
  mutate(taxa = ifelse(taxa %in% bact.high.taxa, taxa, "less than 1% each")) %>%
  group_by(taxa, sample, Sample_type) %>%
  summarize(RA = sum(RA)) %>%
  as.data.frame %>%
  inner_join(new_sample_types)

bact.phyla.high$taxa = factor(bact.phyla.high$taxa, levels=c(sort(bact.high.taxa), "less than 1% each"))
bact.phyla.high$Sample_type_new = factor(bact.phyla.high$Sample_type_new, levels=c("Flowers", "Leaves", "Root tissue", "Rhizosphere soil", "Bulk soil"))

bact.posthoc.res.y = bact.phyla.high %>%
  group_by(taxa) %>%
  summarize(max_RA = max(RA)) %>%
  as.data.frame %>%
  mutate(y = max_RA*1.2) %>%
  inner_join(bact.posthoc.res, by = c("taxa")) %>%
  inner_join(new_sample_types)

bact.posthoc.res.y$taxa = factor(bact.posthoc.res.y$taxa, levels=c(sort(bact.high.taxa), "less than 1% each"))
bact.posthoc.res.y$Sample_type_new = factor(bact.posthoc.res.y$Sample_type_new, levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))
bact.phyla.high$Sample_type_new = factor(bact.phyla.high$Sample_type_new, levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))

bact.comm.plot = ggplot(data=bact.phyla.high, aes(x=Sample_type_new, y=RA, fill=Sample_type_new)) +
  geom_boxplot() +
  geom_text(data=bact.posthoc.res.y, aes(x=Sample_type_new, y=y, label=group)) +
  scale_fill_manual(values=type.col) +
  labs(x="Compartment", y="Relative abundance (% of reads)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  facet_wrap(~taxa, scales = "free_y")
bact.comm.plot

#ggsave(plot = bact.comm.plot, filename = "/home/sam/notebooks/hemp_microbiome/figures/bacteria_taxa_abd.tiff",
#       device = "tiff", width=7, height=7, units="in")

Fungi

fungi.phyla.sum = data.frame()

for (samtype in c("bulk_soil", "rhizosphere_soil", "leaves", "flowers")){
  sub.physeq = subset_samples(physeqITS.fungi.rare, Sample_type == samtype)
  sub.physeq = prune_taxa(taxa_sums(sub.physeq) > 0, sub.physeq)
  sub.otu.table = otu_table(sub.physeq)
  sub.tax.table = data.frame(tax_table(sub.physeq), stringsAsFactors = FALSE) %>%
    rownames_to_column(var="OTU")
  sub.otu.df = data.frame(sub.otu.table) %>%
    rownames_to_column(var="OTU") %>%
    gather(key="sample", value="count", -OTU) %>% 
    filter(count > 0) %>%
    inner_join(sub.tax.table, by="OTU")
  
  sub.phyla.sum = sub.otu.df %>%
    mutate(taxa = ifelse(is.na(Class), ifelse(is.na(Phylum), "Unclassified Fungi", paste("Unclassified", as.character(Phylum))), as.character(Class))) %>%
    group_by(taxa, sample) %>%
    summarize(RA = sum(count)*100) %>%
    as.data.frame %>%
    tidyr::spread(key=sample, value=RA) %>%
    tidyr::gather(key=sample, value=RA, -taxa) %>%
    mutate(RA = ifelse(is.na(RA), 0, RA)) %>%
    mutate(Sample_type = samtype)
  
  fungi.phyla.sum = rbind(fungi.phyla.sum, sub.phyla.sum)
}

fungi.phyla.sum.mean = fungi.phyla.sum %>%
  group_by(taxa, Sample_type) %>%
  summarize(mean_RA = mean(RA)) %>%
  as.data.frame
  
fungi.high.taxa = unique(fungi.phyla.sum.mean[fungi.phyla.sum.mean$mean_RA >= 1 & !(fungi.phyla.sum.mean$taxa %in% c("Unclassified Ascomycota", "Unclassified Fungi")),]$taxa)
fungi.taxa.colors = c("#882E72", "#1965B0", "#5289C7", "#7BAFDE", "#4EB265", "#CAE0AB", "#F7F056", "#F4A736","#E8601C", "#DC050C", "#72190E","grey")
names(fungi.taxa.colors) = c(fungi.high.taxa, "Unclassified Ascomycota", "Unclassified Fungi", "less than 1% each")

new_sample_types = data.frame(Sample_type = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"),
                              Sample_type_new = c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))

fungi.phyla.sum.high = fungi.phyla.sum %>%
  mutate(taxa = ifelse(taxa %in% names(fungi.taxa.colors), taxa, "less than 1% each")) %>%
  group_by(taxa, sample, Sample_type) %>%
  summarize(RA = sum(RA)) %>%
  as.data.frame %>%
  group_by(taxa, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA),
            n_samples=n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  inner_join(new_sample_types)
fungi.phyla.sum.high$taxa = factor(fungi.phyla.sum.high$taxa, levels=c(sort(fungi.high.taxa), "Unclassified Ascomycota", "Unclassified Fungi", "less than 1% each"))
fungi.phyla.sum.high$Sample_type_new = factor(fungi.phyla.sum.high$Sample_type_new, levels=c("Flowers", "Leaves", "Rhizosphere soil", "Bulk soil"))


fungi.comm.plot = ggplot(data=fungi.phyla.sum.high, aes(x=taxa, y=mean_RA, fill=taxa)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=mean_RA-SE_RA, ymax=mean_RA+SE_RA), width = 0.5, size=0.25) +
  scale_fill_manual(values=fungi.taxa.colors) +
  labs(x="Class", y="Relative abundance (% of reads)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  facet_wrap(~Sample_type_new, nrow = 5)
fungi.comm.plot

Here is the same data but in table format.

fungi.phyla.sum.high.df = fungi.phyla.sum.high %>%
  dplyr::select(taxa, Sample_type, mean_RA) %>%
  spread(key=Sample_type, value=mean_RA)

kable(fungi.phyla.sum.high.df, "html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
taxa bulk_soil flowers leaves rhizosphere_soil
Agaricomycetes 1.4391769 1.7568994 0.2101242 2.4014190
Dothideomycetes 6.3134931 13.6970291 27.9274117 6.4913358
Eurotiomycetes 6.0627768 0.0678631 0.0334288 5.0040933
Exobasidiomycetes 0.6121386 35.7751470 39.6872015 4.9017601
Leotiomycetes 1.6475645 0.1545770 0.0668577 1.4019648
Microbotryomycetes 0.1400104 2.4732318 1.8338109 0.1876109
Mortierellomycetes 2.0448033 0.0037702 0.0023878 2.3297858
Sordariomycetes 46.2457671 3.3743025 0.5945559 39.0469368
Tremellomycetes 2.8816098 14.2964862 17.0678128 4.4924273
Unclassified Ascomycota 1.1461318 22.6059418 5.3963706 1.3849093
Unclassified Fungi 29.4445168 5.1915247 6.0840497 29.5060718
less than 1% each 2.0220109 0.6032273 1.0959885 2.8516851

Using nonparametric Kruskal test with Dunn posthoc test, lets see if the relative abundances of these taxa differ across plant compartments.

fungi.posthoc.res = data.frame()
for (subtaxa in fungi.high.taxa){
  taxa.model = kruskal.test(RA~Sample_type, data=filter(fungi.phyla.sum, taxa == subtaxa))
  if (p.adjust(taxa.model$p.value, n=length(high.taxa), method="bonferroni") < 0.05){
    posthoc.res = FSA::dunnTest(RA~Sample_type, data=filter(fungi.phyla.sum, taxa == subtaxa))
    posthoc.res = posthoc.res$res
    posthoc.cld = rcompanion::cldList(comparison = posthoc.res$Comparison,
                                      p.value    = posthoc.res$P.adj,
                                      threshold  = 0.05)
  } else{
    posthoc.cld = data.frame(Group = NA, Letter = NA)
  }
  sub.fungi.posthoc.res = data.frame(taxa = subtaxa, X2 = taxa.model$statistic,
                                    Pvalue = taxa.model$p.value,
                                    padj = p.adjust(taxa.model$p.value, n=length(high.taxa), method="bonferroni"),
                                    Sample_type = posthoc.cld$Group,
                                    group = posthoc.cld$Letter,
                                    stringsAsFactors = FALSE)
  fungi.posthoc.res = rbind(fungi.posthoc.res, sub.fungi.posthoc.res)
}
fungi.posthoc.res$Sample_type = factor(fungi.posthoc.res$Sample_type, levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"))
#write.table(fungi.posthoc.res, "/home/sam/notebooks/hemp_microbiome/data/fungi_phyla_across_compartments.txt", quote=FALSE, row.names = FALSE)

kable(fungi.posthoc.res %>%
        arrange(taxa, Sample_type), "html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
taxa X2 Pvalue padj Sample_type group
Agaricomycetes 31.00921 8e-07 1.1e-05 bulk_soil a
Agaricomycetes 31.00921 8e-07 1.1e-05 rhizosphere_soil a
Agaricomycetes 31.00921 8e-07 1.1e-05 leaves b
Agaricomycetes 31.00921 8e-07 1.1e-05 flowers a
Dothideomycetes 48.13840 0e+00 0.0e+00 bulk_soil a
Dothideomycetes 48.13840 0e+00 0.0e+00 rhizosphere_soil a
Dothideomycetes 48.13840 0e+00 0.0e+00 leaves b
Dothideomycetes 48.13840 0e+00 0.0e+00 flowers a
Eurotiomycetes 73.68029 0e+00 0.0e+00 bulk_soil a
Eurotiomycetes 73.68029 0e+00 0.0e+00 rhizosphere_soil a
Eurotiomycetes 73.68029 0e+00 0.0e+00 leaves b
Eurotiomycetes 73.68029 0e+00 0.0e+00 flowers b
Exobasidiomycetes 70.23297 0e+00 0.0e+00 bulk_soil a
Exobasidiomycetes 70.23297 0e+00 0.0e+00 rhizosphere_soil a
Exobasidiomycetes 70.23297 0e+00 0.0e+00 leaves b
Exobasidiomycetes 70.23297 0e+00 0.0e+00 flowers b
Leotiomycetes 58.74971 0e+00 0.0e+00 bulk_soil a
Leotiomycetes 58.74971 0e+00 0.0e+00 rhizosphere_soil a
Leotiomycetes 58.74971 0e+00 0.0e+00 leaves b
Leotiomycetes 58.74971 0e+00 0.0e+00 flowers b
Microbotryomycetes 58.46442 0e+00 0.0e+00 bulk_soil a
Microbotryomycetes 58.46442 0e+00 0.0e+00 rhizosphere_soil a
Microbotryomycetes 58.46442 0e+00 0.0e+00 leaves b
Microbotryomycetes 58.46442 0e+00 0.0e+00 flowers b
Mortierellomycetes 78.45327 0e+00 0.0e+00 bulk_soil a
Mortierellomycetes 78.45327 0e+00 0.0e+00 rhizosphere_soil a
Mortierellomycetes 78.45327 0e+00 0.0e+00 leaves b
Mortierellomycetes 78.45327 0e+00 0.0e+00 flowers b
Sordariomycetes 70.37733 0e+00 0.0e+00 bulk_soil a
Sordariomycetes 70.37733 0e+00 0.0e+00 rhizosphere_soil a
Sordariomycetes 70.37733 0e+00 0.0e+00 leaves b
Sordariomycetes 70.37733 0e+00 0.0e+00 flowers b
Tremellomycetes 55.38536 0e+00 0.0e+00 bulk_soil a
Tremellomycetes 55.38536 0e+00 0.0e+00 rhizosphere_soil a
Tremellomycetes 55.38536 0e+00 0.0e+00 leaves b
Tremellomycetes 55.38536 0e+00 0.0e+00 flowers b

Now plot the results.

fungi.phyla.high = fungi.phyla.sum %>%
  mutate(taxa = ifelse(taxa %in% fungi.high.taxa, taxa, "less than 1% each")) %>%
  group_by(taxa, sample, Sample_type) %>%
  summarize(RA = sum(RA)) %>%
  as.data.frame %>%
  inner_join(new_sample_types)

fungi.phyla.high$taxa = factor(fungi.phyla.high$taxa, levels=c(sort(fungi.high.taxa), "less than 1% each"))
fungi.phyla.high$Sample_type_new = factor(fungi.phyla.high$Sample_type_new, levels=c("Flowers", "Leaves", "Root tissue", "Rhizosphere soil", "Bulk soil"))

fungi.posthoc.res.y = fungi.phyla.high %>%
  group_by(taxa) %>%
  summarize(max_RA = max(RA)) %>%
  as.data.frame %>%
  mutate(y = max_RA*1.2) %>%
  inner_join(fungi.posthoc.res, by = c("taxa")) %>%
  inner_join(new_sample_types)

fungi.posthoc.res.y$taxa = factor(fungi.posthoc.res.y$taxa, levels=c(sort(fungi.high.taxa), "less than 1% each"))
fungi.posthoc.res.y$Sample_type_new = factor(fungi.posthoc.res.y$Sample_type_new, levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))
fungi.phyla.high$Sample_type_new = factor(fungi.phyla.high$Sample_type_new, levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers"))

fungi.comm.plot = ggplot(data=fungi.phyla.high, aes(x=Sample_type_new, y=RA, fill=Sample_type_new)) +
  geom_boxplot() +
  geom_text(data=fungi.posthoc.res.y, aes(x=Sample_type_new, y=y, label=group)) +
  scale_fill_manual(values=type.col) +
  labs(x="Compartment", y="Relative abundance (% of reads)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none") +
  facet_wrap(~taxa, scales = "free_y")
fungi.comm.plot

#ggsave(plot = fungi.comm.plot, filename = "/home/sam/notebooks/hemp_microbiome/figures/fungi_taxa_abd.tiff",
#       device = "tiff", width=7, height=7, units="in")

Ordination and PERMANOVA

With this analysis I will be statistically comparing the community compositions of the samples. I am testing if the variation in community compositions can be statistically explained by either sample type, field or an interaction between the two. For reference, I am using Bray Curtis dissimilarity for this analysis. I am using this metric mostly to stay consistent because I cannot use UniFrac for the fungi as I do not have a phylogenetic tree for that data.

Bacteria

# Use the function ordinate from phyloseq to get the xy coordinates of an NMDS of the distance matrix.
set.seed(4242)
bact.bray.ord = ordinate(physeq16S.bact.rare, "NMDS", distance="bray", trymax=1000)
bact.bray.ord.df = data.frame(bact.bray.ord$points)
bact.bray.ord.df$sample_ID = rownames(bact.bray.ord.df)

# Add the metadata from the phyloseq to the xy coordinates so that you can do some coloring and other asthetics. Lets also change some of the variable entries to make them look better on a figure (remove "_" and capitalize first letters)
bact.metadata = data.frame(sample_data(physeq16S.bact.rare), stringsAsFactors = FALSE) %>%
  mutate(Field = ifelse(Field == "Blight", "Crittenden North", gsub("_", " ", Field))) %>%
  mutate(Sample_type = factor(gsub("_", " ", paste(toupper(substr(Sample_type, 1, 1)), substr(Sample_type, 2, nchar(Sample_type)), sep="")), 
                              levels=c("Bulk soil", "Rhizosphere soil", "Root tissue", "Leaves", "Flowers")))

bact.bray.ord.df = full_join(bact.bray.ord.df, bact.metadata, by="sample_ID")
# Make a distance matrix using the phyloseq distance function. This is based on the function vegdist from package vegan
bact.bray.dist = phyloseq::distance(physeq16S.bact.rare, "bray")

# Run adonis function with distance ~ Sample_type*field as the formula.
set.seed(4242)
bact.bray.adon = adonis(bact.bray.dist ~ Sample_type*Field, as(sample_data(physeq16S.bact.rare), "data.frame"))
print(bact.bray.adon)
## 
## Call:
## adonis(formula = bact.bray.dist ~ Sample_type * Field, data = as(sample_data(physeq16S.bact.rare),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                    Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Sample_type         4    15.272  3.8179 19.5367 0.32310  0.001 ***
## Field               5     3.795  0.7590  3.8839 0.08029  0.001 ***
## Sample_type:Field  20     7.679  0.3840  1.9648 0.16247  0.001 ***
## Residuals         105    20.519  0.1954         0.43413           
## Total             134    47.265                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now I’ll plot the ordination of this data to visualize this difference.

# Plot the ordination
bacteria.ord.plot = ggplot(data=bact.bray.ord.df, aes(x=MDS1, y=MDS2)) +
  geom_point(aes(color=Sample_type, fill=Sample_type, shape=Field), size=3) +
  scale_color_manual(values=type.col) +
  scale_fill_manual(values=type.col) +
  scale_shape_manual(values=field.shape) +
  labs(x="NMDS 1", y="NMDS 2", fill="Compartment", color="Compartment", shape="Field", 
       title=paste("Bacteria (Stress=", round(bact.bray.ord$stress, digits = 4), ")", sep="")) +
  theme_bw() +
  theme(text = element_text(size=16),
        plot.title = element_text(hjust = 0.5))

bacteria.ord.plot

Now I’ll run the PERMANOVA for each plant compartment individually. In this case, field location is the only explanatory variable.

for (plant_part in c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers")){
  sub.physeq = subset_samples(physeq16S.bact.rare, Sample_type == plant_part)
  sub.physeq = prune_taxa(taxa_sums(sub.physeq) > 0, sub.physeq)
  
  # Make a distance matrix using the phyloseq distance function. This is based on the function vegdist from package vegan
  sub.physeq.dist = phyloseq::distance(sub.physeq, "bray")
  
  # Run adonis function with distance ~ Sample_type*field as the formula.
  set.seed(4242)
  sub.physeq.adon = adonis(sub.physeq.dist ~ Field, as(sample_data(sub.physeq), "data.frame"))
  writeLines(paste("Running", plant_part, "samples now"))
  print(sub.physeq.adon)
  writeLines(paste("\n\nAdjusted pvalue for Field =", p.adjust(sub.physeq.adon$aov.tab$`Pr(>F)`[1], method="bonferroni", n=5)))
  writeLines("\n\n-----\n\n")
}
## Running bulk_soil samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)    
## Field      5    2.4153 0.48306  3.0287 0.4309  0.001 ***
## Residuals 20    3.1899 0.15950         0.5691           
## Total     25    5.6053                 1.0000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.005
## 
## 
## -----
## 
## 
## Running rhizosphere_soil samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      5    2.5197 0.50395  3.3471 0.43205  0.001 ***
## Residuals 22    3.3124 0.15056         0.56795           
## Total     27    5.8321                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.005
## 
## 
## -----
## 
## 
## Running root_tissue samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      5    2.2779 0.45558   2.294 0.35325  0.001 ***
## Residuals 21    4.1704 0.19859         0.64675           
## Total     26    6.4483                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.005
## 
## 
## -----
## 
## 
## Running leaves samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## Field      5    1.5370 0.30741  1.6105 0.26795  0.005 **
## Residuals 22    4.1992 0.19087         0.73205          
## Total     27    5.7362                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.025
## 
## 
## -----
## 
## 
## Running flowers samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      5    2.7244 0.54488  1.9296 0.32542  0.001 ***
## Residuals 20    5.6475 0.28238         0.67458           
## Total     25    8.3719                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.005
## 
## 
## -----

Fungi

# Use the function ordinate from phyloseq to get the xy coordinates of an NMDS of the distance matrix.
set.seed(4242)
fungi.bray.ord = ordinate(physeqITS.fungi.rare, "NMDS", distance="bray", trymax=1000)
fungi.bray.ord.df = data.frame(fungi.bray.ord$points)
fungi.bray.ord.df$sample_ID = rownames(fungi.bray.ord.df)

# Add the metadata from the phyloseq to the xy coordinates so that you can do some coloring and other asthetics. Lets also change some of the variable entries to make them look better on a figure (remove "_" and capitalize first letters)
fungi.metadata = data.frame(sample_data(physeqITS.fungi.rare), stringsAsFactors = FALSE) %>%
  mutate(Field = ifelse(Field == "Blight", "Crittenden North", gsub("_", " ", Field))) %>%
  mutate(Sample_type = gsub("_", " ", paste(toupper(substr(Sample_type, 1, 1)), substr(Sample_type, 2, nchar(Sample_type)), sep="")))
fungi.bray.ord.df = full_join(fungi.bray.ord.df, fungi.metadata, by="sample_ID")
# Make a distance matrix using the phyloseq distance function. This is based on the function vegdist from package vegan
fungi.bray.dist = phyloseq::distance(physeqITS.fungi.rare, "bray")

# Run adonis function with distance ~ Sample_type*field as the formula.
set.seed(4242)
fungi.bray.adon = adonis(fungi.bray.dist ~ Sample_type*Field, as(sample_data(physeqITS.fungi.rare), "data.frame"))
print(fungi.bray.adon)
## 
## Call:
## adonis(formula = fungi.bray.dist ~ Sample_type * Field, data = as(sample_data(physeqITS.fungi.rare),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                   Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Sample_type        3   13.1615  4.3872  49.663 0.47902  0.001 ***
## Field              5    3.9457  0.7891   8.933 0.14361  0.001 ***
## Sample_type:Field 13    4.1848  0.3219   3.644 0.15231  0.001 ***
## Residuals         70    6.1838  0.0883         0.22506           
## Total             91   27.4758                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now I’ll plot the ordination of this data to visualize this difference.

# Plot the ordination
fungal.ord.plot = ggplot(data=fungi.bray.ord.df, aes(x=MDS1, y=MDS2)) +
  geom_point(aes(color=Sample_type, fill=Sample_type, shape=Field), size=3) +
  scale_color_manual(values=type.col) +
  scale_fill_manual(values=type.col) +
  scale_shape_manual(values=field.shape) +
  labs(x="NMDS 1", y="NMDS 2", fill="Compartment", color="Compartment", shape="Field", 
       title=paste("Fungi (Stress=", round(fungi.bray.ord$stress, digits = 4), ")", sep="")) +
  theme_bw() +
  theme(text = element_text(size=16),
        plot.title = element_text(hjust = 0.5))
fungal.ord.plot

Now I’ll run the PERMANOVA for each plant compartment individually. In this case, field location is the only explanatory variable.

for (plant_part in c("bulk_soil", "rhizosphere_soil", "leaves", "flowers")){
  sub.physeq = subset_samples(physeqITS.fungi.rare, Sample_type == plant_part)
  sub.physeq = prune_taxa(taxa_sums(sub.physeq) > 0, sub.physeq)
  
  # Make a distance matrix using the phyloseq distance function. This is based on the function vegdist from package vegan
  sub.physeq.dist = phyloseq::distance(sub.physeq, "bray")
  
  # Run adonis function with distance ~ Sample_type*field as the formula.
  set.seed(4242)
  sub.physeq.adon = adonis(sub.physeq.dist ~ Field, as(sample_data(sub.physeq), "data.frame"))
  writeLines(paste("Running", plant_part, "samples now"))
  print(sub.physeq.adon)
  writeLines(paste("\n\nAdjusted pvalue for Field =", p.adjust(sub.physeq.adon$aov.tab$`Pr(>F)`[1], method="bonferroni", n=4)))
  writeLines("\n\n-----\n\n")
}
## Running bulk_soil samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      4    2.8127 0.70318  6.6302 0.60938  0.001 ***
## Residuals 17    1.8030 0.10606         0.39062           
## Total     21    4.6157                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.004
## 
## 
## -----
## 
## 
## Running rhizosphere_soil samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      5    2.6487 0.52973  4.6837 0.60956  0.001 ***
## Residuals 15    1.6965 0.11310         0.39044           
## Total     20    4.3452                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.004
## 
## 
## -----
## 
## 
## Running leaves samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)    
## Field      5    1.2767 0.25533  4.2841 0.4716  0.001 ***
## Residuals 24    1.4304 0.05960         0.5284           
## Total     29    2.7071                 1.0000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.004
## 
## 
## -----
## 
## 
## Running flowers samples now
## 
## Call:
## adonis(formula = sub.physeq.dist ~ Field, data = as(sample_data(sub.physeq),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Field      4    1.3925 0.34812  3.8868 0.52618  0.001 ***
## Residuals 14    1.2539 0.08956         0.47382           
## Total     18    2.6464                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Adjusted pvalue for Field = 0.004
## 
## 
## -----

Ordination plots together

For the publication, we probably want these ordination to be in the same figure.

plottheme = theme(title = element_blank(),
                  axis.text = element_text(size=10),
                  axis.title = element_blank(),
                  legend.text = element_text(size=10),
                  legend.title = element_text(size=12))

# Get legend from bacterial plot
ord.legend = g_legend(bacteria.ord.plot + plottheme)

# Use cowplot::plot_grid to plot both together. Remove the legend from the plots and add in the full legend
ord.plot = cowplot::plot_grid(bacteria.ord.plot + plottheme + theme(legend.position = "none", plot.title = element_blank()), 
                              fungal.ord.plot + ylim(-1.3, 1.2) + plottheme + theme(legend.position = "none", , plot.title = element_blank()), 
                              ord.legend, ncol=3, rel_widths = c(1,1,.6), labels = c("a", "b", ""))

y.grob <- grid::textGrob("NMDS 2", gp=grid::gpar(fontsize=12), rot=90)

x.grob <- grid::textGrob("NMDS 1", gp=grid::gpar(fontsize=12))

final_ord.plot = gridExtra::arrangeGrob(ord.plot, left = y.grob, bottom = x.grob)

gridExtra::grid.arrange(final_ord.plot)

#ggsave(plot = final_ord.plot, filename = "/home/sam/notebooks/hemp_microbiome/figures/ordinations.tiff",
#       device = "tiff", width=7, height=4, units="in")

Part 2: OTUs of interest

There are a couple ways we can identify OTUs of interest and define a core microbiome.

For OTUs of interest we can identify the most abundant OTUs across all plants and all fields. As they are in high abundance, it may mean that they are adapted to survive with hemp or plants in general, are being selected for, or are pathogens. It is important to keep in mind that even if the OTU is in high abundance, it may not mean it is doing anything with the plant. It may just be super high abundant in the environment or at least in a certain location. High abundance of a particular OTU could also just be an artifact of the compositional nature of sequencing. In other words, it might have a high relative abundance but there may just be very few microbes there in the first place.

Bacteria

bact.physeq = transform_sample_counts(physeq16S.bact.rare, function(x) x/sum(x))
bact.meta = data.frame(sample_data(bact.physeq))
bact.tax.table = data.frame(gsub("D_.__", "", tax_table(bact.physeq)), stringsAsFactors = FALSE) %>%
  rownames_to_column(var="OTU")
bact.OTU.rare.sum = data.frame(otu_table(bact.physeq)) %>%
  rownames_to_column(var="OTU") %>%
  gather(key="sample_ID", value="RA", -OTU) %>%
  mutate(sample_ID = gsub("X", "", gsub("\\.", "-", sample_ID))) %>%
  left_join(bact.meta, by="sample_ID") %>%
  group_by(OTU, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA),
            min_RA = min(RA),
            max_RA = max(RA),
            n_samples = n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  unique() %>%
  inner_join(bact.tax.table, by="OTU")

bact.OTU.rare.soil = bact.OTU.rare.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  dplyr::select(OTU, mean_RA, sd_RA, min_RA, max_RA, n_samples, SE_RA) %>%
  dplyr::rename(soil_mean_RA = mean_RA, soil_sd_RA = sd_RA, soil_min_RA = min_RA, soil_max_RA = max_RA, soil_n_samples = n_samples, soil_SE_RA = SE_RA)

Top 5 most abundant OTUS

Here are the top 5 most abundant bacterial OTUs across all plants:

bact.OTU.rare.abd = bact.OTU.rare.sum %>%
  group_by(Sample_type) %>%
  top_n(n=5, wt=mean_RA) %>%
  as.data.frame %>%
  arrange(Sample_type, -mean_RA) %>%
  left_join(bact.OTU.rare.soil, by="OTU")

kable(bact.OTU.rare.abd %>%
        mutate(FC_abd = mean_RA/soil_mean_RA, mean_RA = mean_RA*100, SE_RA = SE_RA*100) %>%
        mutate(mean_RA = round(mean_RA, 2), SE_RA = round(SE_RA, 2), FC_abd = round(FC_abd, 0)) %>%
        mutate(compartment = paste(Sample_type, " (", n_samples, ")", sep=""),
               RA_SE = paste(mean_RA, " (", SE_RA, ")", sep="")) %>%
        dplyr::select(compartment, OTU, RA_SE, FC_abd, Rank1, Rank2, Rank3, Rank4, Rank5, Rank6, Rank7), "html") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
compartment OTU RA_SE FC_abd Rank1 Rank2 Rank3 Rank4 Rank5 Rank6 Rank7
bulk_soil (26) OTU.10 1.77 (0.25) 1 Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Pseudarthrobacter NA
bulk_soil (26) OTU.28 1.67 (0.14) 1 Bacteria Chloroflexi KD4-96 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
bulk_soil (26) OTU.21 1.29 (0.13) 1 Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Janibacter NA
bulk_soil (26) OTU.66 1.23 (0.23) 1 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
bulk_soil (26) OTU.39 1.12 (0.14) 1 Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Lysobacter NA
flowers (26) OTU.7 12.56 (3.19) 58 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
flowers (26) OTU.4 11.67 (4.41) 2453 Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Lactococcus Lactococcus lactis
flowers (26) OTU.12 6.33 (2.01) 285 Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Pantoea Ambiguous_taxa
flowers (26) OTU.20 6.23 (2) 80 Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Enterobacter NA
flowers (26) OTU.6 5.53 (3.59) 24 Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Bacillus cereus
leaves (28) OTU.8 13.14 (1.88) 237 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
leaves (28) OTU.5 12.81 (1.99) 278 Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
leaves (28) OTU.9 8.55 (1.73) 154 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa
leaves (28) OTU.11 6.53 (1.09) 71 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves (28) OTU.23 5.63 (1.78) 355 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
rhizosphere_soil (28) OTU.10 2 (0.22) 1 Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Pseudarthrobacter NA
rhizosphere_soil (28) OTU.28 1.63 (0.12) 1 Bacteria Chloroflexi KD4-96 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
rhizosphere_soil (28) OTU.21 1.2 (0.16) 1 Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Janibacter NA
rhizosphere_soil (28) OTU.66 1 (0.11) 1 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
rhizosphere_soil (28) OTU.63 0.85 (0.07) 1 Bacteria Acidobacteria Subgroup 6 NA NA NA NA
root_tissue (27) OTU.24 5.69 (0.79) 37 Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces Ambiguous_taxa
root_tissue (27) OTU.17 4.94 (1.33) 1038 Bacteria Proteobacteria Gammaproteobacteria Cellvibrionales Halieaceae NA NA
root_tissue (27) OTU.49 3.47 (0.65) 13 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae NA NA
root_tissue (27) OTU.18 2.68 (0.6) 8 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
root_tissue (27) OTU.117 2.43 (0.45) 14 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Aquabacterium NA

Now I’ll plot this data.

bact.OTU.abd4fig = rbind(bact.OTU.rare.abd %>%
                           filter(Sample_type != "bulk_soil") %>%
                           mutate(Sample_type = stringr::str_to_sentence(gsub("_", " ", Sample_type)),
                                  mean_RA = 100*mean_RA,
                                  SE_RA = 100*SE_RA) %>%
                           dplyr::select(OTU, Sample_type, mean_RA, SE_RA) %>%
                           mutate(part = "Plant"),
                         bact.OTU.rare.abd %>%
                           filter(Sample_type != "bulk_soil") %>%
                           mutate(Sample_type = stringr::str_to_sentence(gsub("_", " ", Sample_type)),
                                  soil_mean_RA = 100*soil_mean_RA,
                                  soil_SE_RA = 100*soil_SE_RA) %>%
                           dplyr::select(OTU, Sample_type, soil_mean_RA, soil_SE_RA) %>%
                           dplyr::rename(mean_RA = soil_mean_RA, SE_RA = soil_SE_RA) %>%
                           mutate(part = "Bulk Soil")) %>%
  mutate(part = factor(part, levels=c("Plant", "Bulk Soil")),
         Sample_type = factor(Sample_type, levels=c("Rhizosphere soil", "Root tissue", "Leaves", "Flowers")))
                        
top_bact.plot = ggplot(data=bact.OTU.abd4fig, aes(x=reorder_within(OTU, -mean_RA, Sample_type), fill=part)) +
  geom_bar(aes(y=mean_RA), stat="identity", color="black", position = "dodge") +
  geom_errorbar(aes(ymin=mean_RA-SE_RA, ymax=mean_RA+SE_RA), width=0.4, position = position_dodge(width = 0.9)) +
  labs(x="Top 5 most abundant bacterial OTUs", y="Relative abundance (%)") +
  scale_fill_manual(values=c("Plant" = "grey", "Bulk soil"="white")) +
  scale_x_reordered() +
  theme_bw() +
  theme(axis.text.x=element_text(size=10, angle=45, hjust=1),
        axis.text.y=element_text(size=10),
        axis.title=element_text(size=12),
        strip.text=element_text(size=12),
        legend.text=element_text(size=10),
        legend.title=element_blank(),
        legend.position = "top") +
  facet_grid(~Sample_type, scales="free_x")

top_bact.plot

Ubiquitous OTUs

The bacteria found on all plants are:

bact.physeq = transform_sample_counts(physeq16S.bact, function(x) x/sum(x))
bact.meta = data.frame(sample_data(bact.physeq))
bact.tax.table = data.frame(gsub("D_.__", "", tax_table(bact.physeq)), stringsAsFactors = FALSE) %>%
  rownames_to_column(var="OTU")
bact.OTU.sum = data.frame(otu_table(bact.physeq)) %>%
  rownames_to_column(var="OTU") %>%
  gather(key="sample_ID", value="RA", -OTU) %>%
  mutate(sample_ID = gsub("X", "", gsub("\\.", "-", sample_ID))) %>%
  left_join(bact.meta, by="sample_ID") %>%
  group_by(OTU, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA),
            min_RA = min(RA),
            max_RA = max(RA),
            n_samples = n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  unique() %>%
  inner_join(bact.tax.table, by="OTU")

bact.OTU.soil = bact.OTU.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  dplyr::select(OTU, mean_RA, sd_RA, min_RA, max_RA, n_samples, SE_RA) %>%
  dplyr::rename(soil_mean_RA = mean_RA, soil_sd_RA = sd_RA, soil_min_RA = min_RA, soil_max_RA = max_RA, soil_n_samples = n_samples, soil_SE_RA = SE_RA)

bact.OTU.ubiq = bact.OTU.sum %>%
  filter(min_RA > 0) %>%
  arrange(Sample_type, -mean_RA)

soil.ubiq = unique(bact.OTU.ubiq[bact.OTU.ubiq$Sample_type == "bulk_soil",]$OTU)

bact.OTU.ubiq = bact.OTU.ubiq %>%
  filter(Sample_type != "bulk_soil") %>%
  mutate(Ubiq_in_soil = ifelse(OTU %in% soil.ubiq, "YES", "NO")) %>%
  left_join(bact.OTU.soil, by="OTU")

kable(bact.OTU.ubiq %>%
        mutate(FC_abd = mean_RA/soil_mean_RA, mean_RA = mean_RA*100, SE_RA = SE_RA*100) %>%
        mutate(mean_RA = round(mean_RA, 2), SE_RA = round(SE_RA, 2), FC_abd = round(FC_abd, 0)) %>%
        mutate(compartment = paste(Sample_type, " (", n_samples, ")", sep=""),
               RA_SE = paste(mean_RA, " (", SE_RA, ")", sep=""),
               Ubiq_in_soil = ifelse(Ubiq_in_soil == "YES", "Y", ifelse(Ubiq_in_soil == "NO", "N", NA))) %>%
        dplyr::select(compartment, OTU, RA_SE, FC_abd, Ubiq_in_soil, Rank1, Rank2, Rank3, Rank4, Rank5, Rank6, Rank7)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
compartment OTU RA_SE FC_abd Ubiq_in_soil Rank1 Rank2 Rank3 Rank4 Rank5 Rank6 Rank7
flowers (26) OTU.6 5.44 (3.58) 24 Y Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Bacillus cereus
flowers (26) OTU.3 5.3 (2.85) 518 N Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
flowers (26) OTU.13 2.89 (0.85) 8 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Burkholderiaceae Ralstonia NA
flowers (26) OTU.14 1.95 (0.46) 4 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
flowers (26) OTU.61 1.22 (0.54) 7 Y Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Ambiguous_taxa
flowers (26) OTU.5 1.01 (0.33) 21 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
flowers (26) OTU.8 0.98 (0.33) 14 Y Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
flowers (26) OTU.9568 0.73 (0.19) 52 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
flowers (26) OTU.66 0.44 (0.15) 0 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
flowers (26) OTU.9 0.38 (0.13) 6 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa
flowers (26) OTU.241 0.29 (0.06) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia Ambiguous_taxa
flowers (26) OTU.37 0.28 (0.08) 1 Y Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
flowers (26) OTU.21 0.24 (0.04) 0 Y Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Janibacter NA
leaves (28) OTU.8 13.05 (1.89) 179 Y Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
leaves (28) OTU.5 12.69 (1.98) 269 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
leaves (28) OTU.9 8.59 (1.7) 142 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa
leaves (28) OTU.11 6.5 (1.08) 79 N Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves (28) OTU.23 5.63 (1.74) 464 N Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves (28) OTU.16 4.89 (0.8) 64 Y Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
leaves (28) OTU.37 2.31 (0.37) 10 Y Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
leaves (28) OTU.25 1.77 (0.48) 124 N Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves (28) OTU.3 1.61 (0.37) 158 N Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
leaves (28) OTU.13 1.26 (0.49) 4 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Burkholderiaceae Ralstonia NA
leaves (28) OTU.18 0.75 (0.26) 2 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
leaves (28) OTU.14 0.72 (0.24) 1 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
leaves (28) OTU.137 0.52 (0.1) 322 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
leaves (28) OTU.20 0.51 (0.21) 5 N Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Enterobacter NA
leaves (28) OTU.52 0.43 (0.14) 1 N Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Nocardioides NA
leaves (28) OTU.1896 0.22 (0.05) 70 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA
leaves (28) OTU.66 0.18 (0.03) 0 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
leaves (28) OTU.28 0.16 (0.05) 0 N Bacteria Chloroflexi KD4-96 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
leaves (28) OTU.6188 0.15 (0.04) 4 N Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae Microbacterium NA
rhizosphere_soil (28) OTU.10 1.98 (0.22) 1 Y Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Pseudarthrobacter NA
rhizosphere_soil (28) OTU.66 1 (0.08) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
rhizosphere_soil (28) OTU.241 0.74 (0.24) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia Ambiguous_taxa
rhizosphere_soil (28) OTU.98 0.71 (0.1) 1 Y Bacteria Acidobacteria Subgroup 6 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
rhizosphere_soil (28) OTU.14 0.55 (0.05) 1 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
rhizosphere_soil (28) OTU.38 0.49 (0.07) 1 N Bacteria Verrucomicrobia Spartobacteria Chthoniobacterales Xiphinematobacteraceae Candidatus Xiphinematobacter uncultured bacterium
rhizosphere_soil (28) OTU.103 0.49 (0.05) 1 N Bacteria Acidobacteria Subgroup 6 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
rhizosphere_soil (28) OTU.18 0.47 (0.08) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
rhizosphere_soil (28) OTU.19 0.41 (0.03) 1 N Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA
rhizosphere_soil (28) OTU.6 0.34 (0.17) 2 Y Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Bacillus cereus
rhizosphere_soil (28) OTU.187 0.34 (0.06) 1 Y Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA
rhizosphere_soil (28) OTU.87 0.33 (0.05) 1 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Microvirga NA
rhizosphere_soil (28) OTU.379 0.24 (0.08) 1 N Bacteria Acidobacteria Subgroup 6 NA NA NA NA
rhizosphere_soil (28) OTU.163 0.22 (0.03) 1 N Bacteria Actinobacteria Thermoleophilia Solirubrobacterales Elev-16S-1332 uncultured bacterium uncultured bacterium
rhizosphere_soil (28) OTU.37 0.2 (0.03) 1 Y Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
rhizosphere_soil (28) OTU.281 0.19 (0.08) 1 Y Bacteria Firmicutes Bacilli Bacillales Planococcaceae NA NA
rhizosphere_soil (28) OTU.437 0.17 (0.02) 2 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Xanthobacteraceae Variibacter NA
rhizosphere_soil (28) OTU.61 0.17 (0.06) 1 Y Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Ambiguous_taxa
rhizosphere_soil (28) OTU.7 0.16 (0.09) 1 N Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
rhizosphere_soil (28) OTU.36 0.14 (0.03) 1 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Devosia NA
rhizosphere_soil (28) OTU.9031 0.13 (0.06) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia uncultured bacterium
rhizosphere_soil (28) OTU.7089 0.11 (0.02) 1 N Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Tetrasphaera uncultured bacterium
rhizosphere_soil (28) OTU.395 0.09 (0.02) 1 N Bacteria Actinobacteria Thermoleophilia Gaiellales uncultured uncultured bacterium uncultured bacterium
rhizosphere_soil (28) OTU.74 0.09 (0.02) 1 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium NA
rhizosphere_soil (28) OTU.8 0.08 (0.03) 1 Y Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
rhizosphere_soil (28) OTU.1408 0.06 (0.01) 1 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
rhizosphere_soil (28) OTU.6032 0.06 (0.01) 1 N Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
root_tissue (27) OTU.49 3.47 (0.64) 14 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae NA NA
root_tissue (27) OTU.18 2.6 (0.58) 8 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
root_tissue (27) OTU.117 2.38 (0.41) 14 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Aquabacterium NA
root_tissue (27) OTU.91 1.11 (0.37) 12 N Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
root_tissue (27) OTU.74 1.1 (0.23) 8 N Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium NA
root_tissue (27) OTU.241 1.05 (0.26) 2 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia Ambiguous_taxa
root_tissue (27) OTU.66 0.93 (0.13) 1 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
root_tissue (27) OTU.14 0.86 (0.11) 2 Y Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
root_tissue (27) OTU.8 0.61 (0.11) 8 Y Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
root_tissue (27) OTU.13 0.57 (0.23) 2 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Burkholderiaceae Ralstonia NA
root_tissue (27) OTU.61 0.32 (0.09) 2 Y Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Ambiguous_taxa
root_tissue (27) OTU.37 0.3 (0.07) 1 Y Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
root_tissue (27) OTU.9031 0.14 (0.04) 2 Y Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia uncultured bacterium
root_tissue (27) OTU.15004 0.14 (0.06) 2 N Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae Microbacterium Ambiguous_taxa

Fungi

fungi.physeq = transform_sample_counts(physeqITS.fungi.rare, function(x) x/sum(x))
fungi.meta = data.frame(sample_data(fungi.physeq))
fungi.tax.table = data.frame(tax_table(fungi.physeq), stringsAsFactors = FALSE) %>%
  rownames_to_column(var="OTU")

fungi.OTU.rare.sum = data.frame(otu_table(fungi.physeq)) %>%
  rownames_to_column(var="OTU") %>%
  gather(key="sample_ID", value="RA", -OTU) %>%
  mutate(sample_ID = gsub("X", "", gsub("\\.", "-", sample_ID))) %>%
  left_join(fungi.meta, by="sample_ID") %>%
  group_by(OTU, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA),
            min_RA = min(RA),
            max_RA = max(RA),
            n_samples = n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  unique() %>%
  inner_join(fungi.tax.table, by="OTU")

fungi.OTU.rare.soil = fungi.OTU.rare.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  dplyr::select(OTU, mean_RA, sd_RA, min_RA, max_RA, n_samples, SE_RA) %>%
  dplyr::rename(soil_mean_RA = mean_RA, soil_sd_RA = sd_RA, soil_min_RA = min_RA, soil_max_RA = max_RA, soil_n_samples = n_samples, soil_SE_RA = SE_RA)

Top 5 most abundant OTUs

Here are the top 5 most abundant fungal OTUs across all plants:

fungi.OTU.rare.abd = fungi.OTU.rare.sum %>%
  group_by(Sample_type) %>%
  top_n(n=5, wt=mean_RA) %>%
  as.data.frame %>%
  arrange(Sample_type, -mean_RA) %>%
  left_join(fungi.OTU.rare.soil, by="OTU")

kable(fungi.OTU.rare.abd %>%
        mutate(FC_abd = mean_RA/soil_mean_RA, mean_RA = mean_RA*100, SE_RA = SE_RA*100) %>%
        mutate(mean_RA = round(mean_RA, 2), SE_RA = round(SE_RA, 2), FC_abd = round(FC_abd, 0)) %>%
        mutate(compartment = paste(Sample_type, " (", n_samples, ")", sep=""),
               RA_SE = paste(mean_RA, " (", SE_RA, ")", sep="")) %>%
        dplyr::select(compartment, OTU, RA_SE, FC_abd, Kingdom, Phylum, Class, Order, Family, Genus, Species)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
compartment OTU RA_SE FC_abd Kingdom Phylum Class Order Family Genus Species
bulk_soil (22) OTU_13 13.97 (2.28) 1 Fungi Ascomycota Sordariomycetes Glomerellales Plectosphaerellaceae Verticillium Verticillium_dahliae_SH466940.07FU
bulk_soil (22) OTU_8 6.93 (0.98) 1 Fungi Ascomycota Sordariomycetes Hypocreales Nectriaceae Fusarium NA
bulk_soil (22) OTU_34 6.84 (1.5) 1 Fungi NA NA NA NA NA NA
bulk_soil (22) OTU_309 5.08 (1.52) 1 Fungi NA NA NA NA NA NA
bulk_soil (22) OTU_33 4.79 (2.07) 1 Fungi Ascomycota Sordariomycetes Hypocreales Hypocreales_fam_Incertae_sedis NA NA
flowers (19) OTU_21 32.3 (5.16) 53 Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
flowers (19) OTU_115 22.33 (3.35) 35 Fungi Ascomycota NA NA NA NA NA
flowers (19) OTU_27 10.8 (2.24) 128 Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
flowers (19) OTU_174 8.18 (1.99) 39 Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
flowers (19) OTU_3520 3.22 (2.26) Inf Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
leaves (30) OTU_21 36.25 (3.11) 60 Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
leaves (30) OTU_174 16.75 (1.67) 80 Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
leaves (30) OTU_27 11.64 (1.64) 138 Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
leaves (30) OTU_115 5.27 (0.68) 8 Fungi Ascomycota NA NA NA NA NA
leaves (30) OTU_125 4.15 (0.59) 12 Fungi NA NA NA NA NA NA
rhizosphere_soil (21) OTU_13 10.94 (1.6) 1 Fungi Ascomycota Sordariomycetes Glomerellales Plectosphaerellaceae Verticillium Verticillium_dahliae_SH466940.07FU
rhizosphere_soil (21) OTU_34 9.57 (1.81) 1 Fungi NA NA NA NA NA NA
rhizosphere_soil (21) OTU_8 7.63 (1.57) 1 Fungi Ascomycota Sordariomycetes Hypocreales Nectriaceae Fusarium NA
rhizosphere_soil (21) OTU_21 4.66 (1.45) 8 Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
rhizosphere_soil (21) OTU_38 4.31 (0.69) 1 Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Didymella Didymella_dimorpha_SH245093.07FU

Now I’ll plot this data.

fungi.OTU.abd4fig = rbind(fungi.OTU.rare.abd %>%
                            filter(Sample_type != "bulk_soil") %>%
                            mutate(Sample_type = stringr::str_to_sentence(gsub("_", " ", Sample_type)),
                                   mean_RA = 100*mean_RA,
                                   SE_RA = 100*SE_RA) %>%
                            dplyr::select(OTU, Sample_type, mean_RA, SE_RA) %>%
                            mutate(part = "Plant"),
                          fungi.OTU.rare.abd %>%
                            filter(Sample_type != "bulk_soil") %>%
                            mutate(Sample_type = stringr::str_to_sentence(gsub("_", " ", Sample_type)),
                                   soil_mean_RA = 100*soil_mean_RA,
                                   soil_SE_RA = 100*soil_SE_RA) %>%
                            dplyr::select(OTU, Sample_type, soil_mean_RA, soil_SE_RA) %>%
                            dplyr::rename(mean_RA = soil_mean_RA, SE_RA = soil_SE_RA) %>%
                            mutate(part = "Bulk Soil")) %>%
  mutate(part = factor(part, levels=c("Plant", "Bulk Soil")),
         Sample_type = factor(Sample_type, levels=c("Rhizosphere soil", "Root tissue", "Leaves", "Flowers")))
                        
top_fungi.plot = ggplot(data=fungi.OTU.abd4fig, aes(x=reorder_within(OTU, -mean_RA, Sample_type), fill=part)) +
  geom_bar(aes(y=mean_RA), stat="identity", color="black", position = "dodge") +
  geom_errorbar(aes(ymin=mean_RA-SE_RA, ymax=mean_RA+SE_RA), width=0.4, position = position_dodge(width = 0.9)) +
  labs(x="Top 5 most abundant fungal OTUs", y="Relative abundance (%)") +
  scale_fill_manual(values=c("Plant" = "grey", "Bulk soil"="white")) +
  scale_x_reordered() +
  theme_bw() +
  theme(axis.text.x=element_text(size=10, angle=45, hjust=1),
        axis.text.y=element_text(size=10),
        axis.title=element_text(size=12),
        strip.text=element_text(size=12),
        legend.text=element_text(size=10),
        legend.title=element_blank(),
        legend.position = "top") +
  facet_grid(~Sample_type, scales="free_x")

top_fungi.plot

Ubiquitous OTUs

The fungi found on all plants are:

fungi.physeq = transform_sample_counts(physeqITS.fungi, function(x) x/sum(x))
fungi.meta = data.frame(sample_data(fungi.physeq))
fungi.tax.table = data.frame(tax_table(fungi.physeq), stringsAsFactors = FALSE) %>%
  rownames_to_column(var="OTU")

fungi.OTU.sum = data.frame(otu_table(fungi.physeq)) %>%
  rownames_to_column(var="OTU") %>%
  gather(key="sample_ID", value="RA", -OTU) %>%
  mutate(sample_ID = gsub("X", "", gsub("\\.", "-", sample_ID))) %>%
  left_join(fungi.meta, by="sample_ID") %>%
  group_by(OTU, Sample_type) %>%
  summarize(mean_RA = mean(RA),
            sd_RA = sd(RA),
            min_RA = min(RA),
            max_RA = max(RA),
            n_samples = n()) %>%
  as.data.frame %>%
  mutate(SE_RA = sd_RA/sqrt(n_samples)) %>%
  unique() %>%
  inner_join(fungi.tax.table, by="OTU")

fungi.OTU.soil = fungi.OTU.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  dplyr::select(OTU, mean_RA, sd_RA, min_RA, max_RA, n_samples, SE_RA) %>%
  dplyr::rename(soil_mean_RA = mean_RA, soil_sd_RA = sd_RA, soil_min_RA = min_RA, soil_max_RA = max_RA, soil_n_samples = n_samples, soil_SE_RA = SE_RA)

fungi.OTU.ubiq = fungi.OTU.sum %>%
  filter(min_RA > 0) %>%
  arrange(Sample_type, -mean_RA)

soil.ubiq = unique(fungi.OTU.ubiq[fungi.OTU.ubiq$Sample_type == "bulk_soil",]$OTU)

fungi.OTU.ubiq = fungi.OTU.ubiq %>%
  filter(Sample_type != "bulk_soil") %>%
  mutate(Ubiq_in_soil = ifelse(OTU %in% soil.ubiq, "YES", "NO")) %>%
  left_join(fungi.OTU.soil, by="OTU")

kable(fungi.OTU.ubiq %>%
        mutate(FC_abd = mean_RA/soil_mean_RA, mean_RA = mean_RA*100, SE_RA = SE_RA*100) %>%
        mutate(mean_RA = round(mean_RA, 2), SE_RA = round(SE_RA, 2), FC_abd = round(FC_abd, 0)) %>%
        mutate(compartment = paste(Sample_type, " (", n_samples, ")", sep=""), 
               RA_SE = paste(mean_RA, " (", SE_RA, ")", sep=""),
               Ubiq_in_soil = ifelse(Ubiq_in_soil == "YES", "Y", ifelse(Ubiq_in_soil == "NO", "N", NA))) %>%
        dplyr::select(compartment, OTU, RA_SE, FC_abd, Ubiq_in_soil, Kingdom, Phylum, Class, Order, Family, Genus, Species)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
compartment OTU RA_SE FC_abd Ubiq_in_soil Kingdom Phylum Class Order Family Genus Species
flowers (19) OTU_21 31.76 (5.07) 50 N Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
flowers (19) OTU_115 22.66 (3.42) 37 Y Fungi Ascomycota NA NA NA NA NA
flowers (19) OTU_27 10.76 (2.26) 133 N Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
flowers (19) OTU_174 8.23 (2.03) 49 N Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
flowers (19) OTU_125 2.13 (0.45) 6 N Fungi NA NA NA NA NA NA
flowers (19) OTU_333 1.47 (0.3) 269 N Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU
flowers (19) OTU_133 0.73 (0.27) 14 N Fungi Basidiomycota Agaricomycetes Russulales NA NA NA
flowers (19) OTU_360 0.25 (0.04) 272 N Fungi Ascomycota Dothideomycetes Capnodiales Mycosphaerellaceae NA NA
leaves (30) OTU_21 35.64 (3.07) 56 N Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
leaves (30) OTU_174 16.98 (1.75) 102 N Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
leaves (30) OTU_27 11.73 (1.59) 145 N Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
leaves (30) OTU_115 5.11 (0.69) 8 Y Fungi Ascomycota NA NA NA NA NA
leaves (30) OTU_125 4.18 (0.54) 12 N Fungi NA NA NA NA NA NA
leaves (30) OTU_38 3.52 (0.75) 1 Y Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Didymella Didymella_dimorpha_SH245093.07FU
leaves (30) OTU_84 2.27 (0.48) 63 N Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Neoascochyta NA
leaves (30) OTU_118 1.74 (0.33) 149 N Fungi Ascomycota Dothideomycetes Pleosporales Pleosporaceae Alternaria Alternaria_infectoria_SH434793.07FU
leaves (30) OTU_331 0.85 (0.11) 168 N Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Dioszegia Dioszegia_hungarica_SH182099.07FU
leaves (30) OTU_338 0.79 (0.11) 162 N Fungi Basidiomycota NA NA NA NA NA
leaves (30) OTU_333 0.77 (0.17) 141 N Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU
leaves (30) OTU_346 0.7 (0.11) 402 N Fungi NA NA NA NA NA NA
leaves (30) OTU_359 0.67 (0.07) 71 N Fungi Ascomycota Dothideomycetes Pleosporales NA NA NA
leaves (30) OTU_332 0.41 (0.09) 22 N Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium NA
leaves (30) OTU_360 0.31 (0.04) 333 N Fungi Ascomycota Dothideomycetes Capnodiales Mycosphaerellaceae NA NA
rhizosphere_soil (21) OTU_13 11.03 (1.61) 1 Y Fungi Ascomycota Sordariomycetes Glomerellales Plectosphaerellaceae Verticillium Verticillium_dahliae_SH466940.07FU
rhizosphere_soil (21) OTU_34 9.46 (1.8) 1 Y Fungi NA NA NA NA NA NA
rhizosphere_soil (21) OTU_8 7.8 (1.5) 1 Y Fungi Ascomycota Sordariomycetes Hypocreales Nectriaceae Fusarium NA
rhizosphere_soil (21) OTU_38 4.21 (0.69) 1 Y Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Didymella Didymella_dimorpha_SH245093.07FU
rhizosphere_soil (21) OTU_5 4.11 (1.05) 1 Y Fungi NA NA NA NA NA NA
rhizosphere_soil (21) OTU_309 2.82 (0.6) 1 Y Fungi NA NA NA NA NA NA
rhizosphere_soil (21) OTU_66 1.55 (0.33) 1 Y Fungi Mortierellomycota Mortierellomycetes Mortierellales Mortierellaceae Mortierella NA
rhizosphere_soil (21) OTU_115 0.92 (0.37) 1 Y Fungi Ascomycota NA NA NA NA NA
rhizosphere_soil (21) OTU_27 0.48 (0.19) 6 N Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU

Plot 5 top abunant OTU figures together

top_abd.leg = g_legend(top_fungi.plot)

top_abd.plot = cowplot::plot_grid(top_abd.leg,
                                  top_bact.plot + theme(legend.position = "none",axis.title = element_blank()),
                                  top_fungi.plot + theme(legend.position = "none",axis.title = element_blank()), 
                                  nrow=3, rel_heights = c(0.2,1,1), labels = c("", "a", "b"))

y.grob <- grid::textGrob("Relative abundance (% of reads)", gp=grid::gpar(fontsize=12), rot=90)

x.grob <- grid::textGrob("Top 5 most abundant OTUs", gp=grid::gpar(fontsize=12))

final_top_abd.plot = gridExtra::arrangeGrob(top_abd.plot, left = y.grob, bottom = x.grob)

gridExtra::grid.arrange(final_top_abd.plot)

#ggsave(plot = final_top_abd.plot, filename = "/home/sam/notebooks/hemp_microbiome/figures/top_abd.tiff",
#       device = "tiff", width=7, height=6.5, units="in")

Part 3: Significantly enriched OTUs

Now lets use DESeq2 to see if any OTUs are enriched in each of the plant compartments compared to the bulk soil. This might indicate that they are interacting in some way with the hemp plants. I will combine these results with those above to identify the core OTUs. Note that for this step I will be using the unrarefied OTU tables.

Bacteria

bact.deseq.df = data.frame()
for (pcomp in c("rhizosphere_soil", "root_tissue", "leaves", "flowers")){
  # Subset the data to just include bulk soil and plant compartment samples.
  physeq16S.bact.BSPC = subset_samples(physeq16S.bact, Sample_type %in% c("bulk_soil", pcomp))
  physeq16S.bact.BSPC = prune_taxa(taxa_sums(physeq16S.bact.BSPC) > 0, physeq16S.bact.BSPC)
  
  bact.metadata = data.frame(sample_data(physeq16S.bact.BSPC), stringsAsFactors = F)
  
  # Set the sample types as a factor
  sample_data(physeq16S.bact.BSPC)$Sample_type = factor(sample_data(physeq16S.bact.BSPC)$Sample_type, levels = c("bulk_soil", pcomp))
  sample_data(physeq16S.bact.BSPC)$Field = as.factor(sample_data(physeq16S.bact.BSPC)$Field)
  
  # Get the taxonomy table
  bact.tax.table = data.frame(gsub("D_.__", "", tax_table(physeq16S.bact.BSPC)), stringsAsFactors = F) %>%
    rownames_to_column(var="OTU") %>%
    dplyr::select(OTU, Rank1, Rank2, Rank3, Rank4, Rank5, Rank6, Rank7)
  colnames(bact.tax.table) = c("OTU", "Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species")
  
  # I'll also remove sparse taxa. In this case, taxa must be found with at least a count of 5 in at least 10 samples.
  OTU.table = as.matrix(otu_table(physeq16S.bact.BSPC))
  OTU.table[OTU.table < 5] = 0
  OTU.table[OTU.table >= 5] = 1
  #sparse = length(sample_names(physeq16S.bact.BSPC)) * 0.25
  sparse = 1
  sampleCount.df = data.frame(sample = rowSums(OTU.table)) %>%
    rownames_to_column(var="OTU") %>%
    filter(sample >= sparse)
  physeq16S.bact.BSPC.sparse = prune_taxa(sampleCount.df$OTU, physeq16S.bact.BSPC)
  
  # Run DESeq2
  bact.deseq = phyloseq_to_deseq2(physeq16S.bact.BSPC.sparse, ~ Field + Sample_type)
  bact.deseq = DESeq(bact.deseq, betaPrior=TRUE)
  bact.deseq.res = results(bact.deseq, lfcThreshold = .25,
                      contrast = c("Sample_type", pcomp, "bulk_soil"), 
                      altHypothesis = "greater",
                      test="Wald")
  bact.deseq.df = rbind(bact.deseq.df, data.frame(bact.deseq.res) %>%
                          rownames_to_column(var="OTU") %>%
                          mutate(Sample_type = pcomp))
}

bact.deseq.df = inner_join(bact.deseq.df, bact.tax.table, by="OTU") %>%
  filter(!(is.na(padj))) %>%
  arrange(-log2FoldChange) %>%
  mutate(sig = ifelse(padj < 0.1, Phylum, "not significant"))
bact.deseq.df$Phylum = factor(bact.deseq.df$Phylum, levels = unique(bact.deseq.df$Phylum))

kable(bact.deseq.df %>% filter(padj < 0.1)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
OTU baseMean log2FoldChange lfcSE stat pvalue padj Sample_type Domain Phylum Class Order Family Genus Species sig
OTU.4 3571.799208 9.920785 0.7234977 13.366712 0.0000000 0.0000000 flowers Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Lactococcus Lactococcus lactis Firmicutes
OTU.17 660.991528 8.623127 0.7208910 11.614969 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Gammaproteobacteria Cellvibrionales Halieaceae NA NA Proteobacteria
OTU.23 634.649652 8.297768 0.6686146 12.036482 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.12 2876.106045 8.034425 0.7291993 10.675306 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Pantoea Ambiguous_taxa Proteobacteria
OTU.5 1420.081052 7.646873 0.4608996 16.048773 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.164 103.511102 7.619215 0.7266164 10.141824 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter Ambiguous_taxa Bacteroidetes
OTU.25 249.704990 7.530690 0.7467881 9.749339 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.9 1289.372673 7.463095 0.5666383 12.729629 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.205 61.213252 7.308270 0.8751016 8.065658 0.0000000 0.0000000 leaves Bacteria Deinococcus-Thermus Deinococci Deinococcales Deinococcaceae Deinococcus Ambiguous_taxa Deinococcus-Thermus
OTU.11 836.188333 7.167521 0.5885840 11.752818 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.259 112.420993 7.165476 0.7844335 8.815886 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Serratia NA Proteobacteria
OTU.11576 90.099056 7.131494 0.7947227 8.658988 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.751 7.562005 6.748306 1.6923285 3.839861 0.0000616 0.0037385 leaves Bacteria Proteobacteria Gammaproteobacteria Orbales Orbaceae Gilliamella NA Proteobacteria
OTU.174 32.420036 6.747010 0.7513153 8.647515 0.0000000 0.0000000 leaves Bacteria Proteobacteria Deltaproteobacteria Myxococcales P3OB-42 uncultured bacterium uncultured bacterium Proteobacteria
OTU.9568 53.295623 6.744988 0.5602651 11.592704 0.0000000 0.0000000 flowers Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA Proteobacteria
OTU.70 137.195911 6.734907 0.7544133 8.595960 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium NA Proteobacteria
OTU.275 41.316155 6.730361 0.7370718 8.792036 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.8 1331.035714 6.694275 0.3908336 16.488539 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium Proteobacteria
OTU.79 253.877489 6.687770 0.5027513 12.805078 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Ambiguous_taxa Ambiguous_taxa Actinobacteria
OTU.3139 17.786374 6.654592 0.9784827 6.545432 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Amycolatopsis NA Actinobacteria
OTU.180 26.950068 6.649146 0.6958561 9.196077 0.0000000 0.0000000 flowers Bacteria Firmicutes Bacilli Bacillales Staphylococcaceae Staphylococcus Ambiguous_taxa Firmicutes
OTU.34 46.623410 6.620527 0.8378273 7.603628 0.0000000 0.0000000 flowers Bacteria Firmicutes Bacilli Bacillales Paenibacillaceae Paenibacillus Ambiguous_taxa Firmicutes
OTU.137 58.961277 6.471813 0.5751450 10.817816 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.7 3381.280443 6.212476 0.5919354 10.072849 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa Proteobacteria
OTU.106 57.794591 6.145799 0.6769505 8.709350 0.0000000 0.0000000 leaves Bacteria Proteobacteria Deltaproteobacteria Myxococcales P3OB-42 uncultured bacterium uncultured bacterium Proteobacteria
OTU.2361 44.630581 6.116103 0.8044843 7.291755 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.4673 181.481231 5.960209 0.7944623 7.187514 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.86 92.489034 5.862703 0.5696949 9.852122 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Caulobacter Ambiguous_taxa Proteobacteria
OTU.20 1546.561799 5.852103 0.6211593 9.018786 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Enterobacter NA Proteobacteria
OTU.266 43.205810 5.737723 0.6224636 8.816135 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales Rhodospirillaceae Dongia Ambiguous_taxa Proteobacteria
OTU.9349 57.635863 5.595116 0.7459712 7.165312 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Pantoea NA Proteobacteria
OTU.143 122.551752 5.586324 0.4876069 10.943906 0.0000000 0.0000000 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Niastella NA Bacteroidetes
OTU.3 201.277688 5.570316 0.6148167 8.653499 0.0000000 0.0000000 leaves Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa Proteobacteria
OTU.357 18.283458 5.520240 0.8361042 6.303329 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium uncultured bacterium Proteobacteria
OTU.194 18.166968 5.461387 0.6837640 7.621617 0.0000000 0.0000000 leaves Bacteria Actinobacteria Actinobacteria Frankiales Geodermatophilaceae Geodermatophilus Ambiguous_taxa Actinobacteria
OTU.6080 51.188522 5.428934 0.6632473 7.808450 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.348 10.055675 5.427031 0.8611814 6.011545 0.0000000 0.0000001 flowers Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Streptococcus uncultured bacterium Firmicutes
OTU.135 70.563954 5.283381 0.6260175 8.040321 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Tahibacter NA Proteobacteria
OTU.8114 16.995512 5.195937 0.7145425 6.921824 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Lentzea NA Actinobacteria
OTU.332 19.446279 5.194073 0.6895378 7.170127 0.0000000 0.0000000 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Cryomorphaceae Fluviicola Ambiguous_taxa Bacteroidetes
OTU.6912 9.480110 5.157315 0.7217388 6.799295 0.0000000 0.0000000 flowers Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae NA NA Proteobacteria
OTU.15206 29.352684 5.133485 0.6474265 7.542917 0.0000000 0.0000000 leaves Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae Rathayibacter NA Actinobacteria
OTU.94 29.364549 5.117718 0.6688757 7.277463 0.0000000 0.0000000 flowers Bacteria Firmicutes Bacilli Bacillales Paenibacillaceae Brevibacillus NA Firmicutes
OTU.4689 13.003557 5.101727 0.8737044 5.553053 0.0000000 0.0000011 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.3696 23.558742 5.055173 0.7557208 6.358398 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.7372 18.627231 5.019477 0.6722214 7.095099 0.0000000 0.0000000 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.16 702.134810 4.982947 0.4407824 10.737604 0.0000000 0.0000000 leaves Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA Actinobacteria
OTU.108 14.340212 4.982034 0.7943707 5.956960 0.0000000 0.0000002 flowers Bacteria Firmicutes Bacilli Bacillales Paenibacillaceae Paenibacillus NA Firmicutes
OTU.231 24.198949 4.972353 1.4898030 3.169784 0.0007628 0.0374194 leaves Bacteria Proteobacteria Gammaproteobacteria Orbales Orbaceae Orbus NA Proteobacteria
OTU.84 120.682549 4.967356 0.4731728 9.969627 0.0000000 0.0000000 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Niastella NA Bacteroidetes
OTU.5 290.678854 4.939744 0.5476515 8.563373 0.0000000 0.0000000 flowers Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.389 12.428940 4.920457 0.6917472 6.751682 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales Acetobacteraceae Roseomonas NA Proteobacteria
OTU.44 173.761221 4.892239 0.5060671 9.173170 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Streptosporangiales Thermomonosporaceae NA NA Actinobacteria
OTU.4541 26.414807 4.823582 0.5988560 7.637199 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Variovorax NA Proteobacteria
OTU.41 114.059842 4.817698 0.4821834 9.472946 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.857 28.153150 4.803418 1.3895366 3.276933 0.0005247 0.0418243 flowers Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Acinetobacter uncultured bacterium Proteobacteria
OTU.93 47.312019 4.796074 0.7761226 5.857417 0.0000000 0.0000002 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Aurantimonadaceae Aureimonas Ambiguous_taxa Proteobacteria
OTU.2320 2.974918 4.772575 1.2363004 3.658152 0.0001270 0.0039395 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales Rhodospirillaceae Ferrovibrio Ambiguous_taxa Proteobacteria
OTU.345 18.953948 4.689072 1.0845472 4.093019 0.0000213 0.0013577 leaves Bacteria Actinobacteria Actinobacteria Kineosporiales Kineosporiaceae NA NA Actinobacteria
OTU.572 7.273330 4.655254 0.8915572 4.941079 0.0000004 0.0000268 leaves Bacteria Proteobacteria Deltaproteobacteria Myxococcales P3OB-42 uncultured bacterium uncultured bacterium Proteobacteria
OTU.237 26.125605 4.636857 0.7521489 5.832432 0.0000000 0.0000002 leaves Bacteria Actinobacteria Actinobacteria Corynebacteriales Nocardiaceae Rhodococcus NA Actinobacteria
OTU.3837 7.433121 4.620387 0.6921786 6.313959 0.0000000 0.0000000 flowers Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Micrococcus NA Actinobacteria
OTU.1896 35.399977 4.602266 0.5299586 8.212465 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA Proteobacteria
OTU.188 14.890436 4.499234 1.3101985 3.243199 0.0005910 0.0456347 flowers Bacteria Firmicutes Bacilli Lactobacillales Enterococcaceae Enterococcus NA Firmicutes
OTU.411 13.210884 4.499042 0.8142147 5.218577 0.0000001 0.0000064 root_tissue Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae uncultured Ambiguous_taxa Bacteroidetes
OTU.650 7.797380 4.498628 1.1022851 3.854383 0.0000580 0.0036094 leaves Bacteria Proteobacteria Deltaproteobacteria Myxococcales P3OB-42 uncultured bacterium uncultured bacterium Proteobacteria
OTU.46 114.231967 4.484545 0.6309144 6.711759 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Amycolatopsis NA Actinobacteria
OTU.48 68.391639 4.465536 0.7093122 5.943132 0.0000000 0.0000001 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter Ambiguous_taxa Bacteroidetes
OTU.99 41.785348 4.456881 1.3152738 3.198483 0.0006908 0.0487679 flowers Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales Acetobacteraceae Asaia NA Proteobacteria
OTU.11 110.981835 4.405180 0.6349852 6.543742 0.0000000 0.0000000 flowers Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.10204 11.044115 4.397589 0.8308867 4.991763 0.0000003 0.0000212 leaves Bacteria Actinobacteria Actinobacteria Corynebacteriales Nocardiaceae Rhodococcus NA Actinobacteria
OTU.5067 8.485922 4.386048 0.6355734 6.507585 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Caulobacter uncultured bacterium Proteobacteria
OTU.670 11.332319 4.301064 0.9729205 4.163818 0.0000156 0.0010236 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Aurantimonadaceae Aureimonas NA Proteobacteria
OTU.330 31.797562 4.276651 0.7065150 5.699315 0.0000000 0.0000005 root_tissue Bacteria Actinobacteria Actinobacteria Micromonosporales Micromonosporaceae Actinoplanes Ambiguous_taxa Actinobacteria
OTU.1030 6.525290 4.243806 0.7684997 5.196887 0.0000001 0.0000069 root_tissue Bacteria Saccharibacteria Ambiguous_taxa Ambiguous_taxa Ambiguous_taxa Ambiguous_taxa Ambiguous_taxa Saccharibacteria
OTU.746 16.502157 4.192488 0.5613301 7.023475 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhodospirillales Rhodospirillaceae Dongia Ambiguous_taxa Proteobacteria
OTU.180 10.720015 4.170523 0.6982414 5.614853 0.0000000 0.0000008 leaves Bacteria Firmicutes Bacilli Bacillales Staphylococcaceae Staphylococcus Ambiguous_taxa Firmicutes
OTU.120 45.321090 4.127426 0.5154983 7.521705 0.0000000 0.0000000 leaves Bacteria Actinobacteria Actinobacteria Kineosporiales Kineosporiaceae Kineococcus Ambiguous_taxa Actinobacteria
OTU.586 7.890956 4.109769 0.8368538 4.612238 0.0000020 0.0000913 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Chitinophaga Ambiguous_taxa Bacteroidetes
OTU.24 650.741983 4.099233 0.3756662 10.246419 0.0000000 0.0000000 root_tissue Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces Ambiguous_taxa Actinobacteria
OTU.269 268.591616 4.096492 0.9599603 4.006928 0.0000308 0.0030400 flowers Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA Proteobacteria
OTU.8 415.309449 4.070808 0.4632500 8.247832 0.0000000 0.0000000 flowers Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium Proteobacteria
OTU.253 20.948096 4.038982 0.6096549 6.214962 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadales Incertae Sedis NA NA Proteobacteria
OTU.23 34.911208 4.015739 0.7741019 4.864656 0.0000006 0.0000708 flowers Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.12779 4.595123 4.008788 0.8130714 4.622949 0.0000019 0.0000887 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Luteibacter NA Proteobacteria
OTU.207 24.781328 3.997731 0.8065566 4.646581 0.0000017 0.0000809 root_tissue Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Dyadobacter NA Bacteroidetes
OTU.3153 2.935691 3.964226 0.7640261 4.861386 0.0000006 0.0000332 root_tissue Bacteria Saccharibacteria uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium Saccharibacteria
OTU.226 32.214272 3.953171 0.5056737 7.323244 0.0000000 0.0000000 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae NA NA Bacteroidetes
OTU.107 31.990477 3.949577 0.8597073 4.303298 0.0000084 0.0003169 root_tissue Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Acinetobacter uncultured bacterium Proteobacteria
OTU.1299 5.456008 3.915910 0.7253159 5.054225 0.0000002 0.0000134 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Taibaiella NA Bacteroidetes
OTU.216 7.664209 3.899453 1.0786581 3.383327 0.0003581 0.0194346 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter NA Bacteroidetes
OTU.5793 21.301971 3.890938 0.8115524 4.486387 0.0000036 0.0001469 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Stenotrophomonas uncultured bacterium Proteobacteria
OTU.2524 34.124376 3.878350 0.7172057 5.059010 0.0000002 0.0000134 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae NA NA Proteobacteria
OTU.197 19.075806 3.845111 0.6892690 5.215832 0.0000001 0.0000064 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Dokdonella NA Proteobacteria
OTU.244 13.270226 3.816975 0.9727477 3.666907 0.0001228 0.0112340 flowers Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Enhydrobacter NA Proteobacteria
OTU.797 4.856738 3.799665 0.9785058 3.627638 0.0001430 0.0126209 flowers Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Streptococcus uncultured bacterium Firmicutes
OTU.50 102.029310 3.720239 0.5265826 6.590113 0.0000000 0.0000000 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium NA Proteobacteria
OTU.629 2.487706 3.715496 0.9831157 3.525014 0.0002117 0.0174395 flowers Bacteria Actinobacteria Actinobacteria Corynebacteriales Corynebacteriaceae Lawsonella NA Actinobacteria
OTU.1243 3.918423 3.711158 1.1611450 2.980815 0.0014374 0.0329512 root_tissue Bacteria Proteobacteria Gammaproteobacteria HTA4 uncultured bacterium uncultured bacterium uncultured bacterium Proteobacteria
OTU.763 7.118758 3.701967 0.6537577 5.280194 0.0000001 0.0000049 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Asticcacaulis Ambiguous_taxa Proteobacteria
OTU.1308 3.517263 3.699751 1.1061419 3.118723 0.0009082 0.0228019 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae uncultured uncultured bacterium Proteobacteria
OTU.7870 3.934522 3.687898 0.9148726 3.757789 0.0000857 0.0050848 leaves Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium Bacteroidetes
OTU.3717 13.322033 3.657838 0.6834608 4.986150 0.0000003 0.0000180 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Niastella Ambiguous_taxa Bacteroidetes
OTU.1493 6.953740 3.631323 0.8384225 4.032958 0.0000275 0.0009368 root_tissue Bacteria Proteobacteria Deltaproteobacteria Myxococcales Sandaracinaceae uncultured uncultured bacterium Proteobacteria
OTU.585 7.372631 3.624998 0.7160717 4.713211 0.0000012 0.0000628 root_tissue Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Chryseolinea Ambiguous_taxa Bacteroidetes
OTU.70 6.711629 3.616253 0.8841161 3.807479 0.0000702 0.0066713 flowers Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium NA Proteobacteria
OTU.237 4.107752 3.600160 0.7863218 4.260546 0.0000102 0.0010954 flowers Bacteria Actinobacteria Actinobacteria Corynebacteriales Nocardiaceae Rhodococcus NA Actinobacteria
OTU.3588 24.848295 3.570163 0.6504389 5.104496 0.0000002 0.0000109 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Saccharothrix NA Actinobacteria
OTU.8419 4.477717 3.567416 0.7038892 4.712981 0.0000012 0.0000628 root_tissue Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae NA NA Actinobacteria
OTU.90 25.626388 3.562114 0.7255389 4.565040 0.0000025 0.0002805 flowers Bacteria Firmicutes Bacilli Bacillales Family XII Exiguobacterium NA Firmicutes
OTU.6513 6.918031 3.533487 0.9130439 3.596198 0.0001615 0.0093605 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.45 196.840824 3.521574 0.4068931 8.040377 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium NA Proteobacteria
OTU.505 3.246543 3.520570 0.9439773 3.464670 0.0002654 0.0077752 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Flavobacteriaceae Chryseobacterium Ambiguous_taxa Bacteroidetes
OTU.89 66.468401 3.495292 0.6148078 5.278547 0.0000001 0.0000049 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Dokdonella NA Proteobacteria
OTU.517 6.752405 3.494807 0.7730227 4.197558 0.0000135 0.0004905 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Chitinophaga NA Bacteroidetes
OTU.254 30.884275 3.488399 0.6918402 4.680847 0.0000014 0.0000704 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Tahibacter uncultured bacterium Proteobacteria
OTU.7395 23.035251 3.482737 0.5689868 5.681568 0.0000000 0.0000006 root_tissue Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae Ambiguous_taxa Ambiguous_taxa Actinobacteria
OTU.204 39.733921 3.479703 0.6439437 5.015505 0.0000003 0.0000159 root_tissue Bacteria Proteobacteria Betaproteobacteria Nitrosomonadales Gallionellaceae NA NA Proteobacteria
OTU.5756 9.838949 3.474829 0.7011108 4.599599 0.0000021 0.0000930 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Caulobacter NA Proteobacteria
OTU.10729 5.026170 3.460692 0.9237574 3.475687 0.0002548 0.0075678 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Aquabacterium Ambiguous_taxa Proteobacteria
OTU.205 2.976369 3.437963 0.9939876 3.207247 0.0006701 0.0486976 flowers Bacteria Deinococcus-Thermus Deinococci Deinococcales Deinococcaceae Deinococcus Ambiguous_taxa Deinococcus-Thermus
OTU.9568 18.935697 3.412282 0.6491354 4.871528 0.0000006 0.0000372 leaves Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA Proteobacteria
OTU.9 41.263112 3.410979 0.5379818 5.875625 0.0000000 0.0000003 flowers Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa Proteobacteria
OTU.1125 3.162377 3.373572 1.0580100 2.952309 0.0015770 0.0974214 flowers Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Rothia uncultured bacterium Actinobacteria
OTU.34 7.280941 3.356052 0.9591352 3.238388 0.0006010 0.0306649 leaves Bacteria Firmicutes Bacilli Bacillales Paenibacillaceae Paenibacillus Ambiguous_taxa Firmicutes
OTU.42 168.413813 3.330700 0.4580325 6.725942 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae NA NA Proteobacteria
OTU.1067 5.287667 3.306502 1.0100764 3.026011 0.0012390 0.0296233 root_tissue Bacteria Armatimonadetes Armatimonadia Armatimonadales uncultured bacterium uncultured bacterium uncultured bacterium Armatimonadetes
OTU.1665 2.507678 3.303193 1.1002090 2.775103 0.0027592 0.0564969 root_tissue Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae uncultured NA Bacteroidetes
OTU.1080 5.027603 3.268440 1.0412091 2.898976 0.0018719 0.0402844 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Cryomorphaceae Fluviicola uncultured bacterium Bacteroidetes
OTU.4689 3.605243 3.259106 0.9269858 3.246119 0.0005850 0.0156160 root_tissue Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.107 16.865028 3.249187 0.8333254 3.599059 0.0001597 0.0136063 flowers Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Acinetobacter uncultured bacterium Proteobacteria
OTU.137 4.671164 3.248178 0.7453162 4.022694 0.0000288 0.0029619 flowers Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.858 4.452337 3.181725 0.9813983 2.987293 0.0014073 0.0326153 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Flavobacteriaceae Chryseobacterium NA Bacteroidetes
OTU.1253 8.076241 3.172270 0.8510477 3.433733 0.0002977 0.0085997 root_tissue Bacteria Proteobacteria Betaproteobacteria Nitrosomonadales Nitrosomonadaceae uncultured uncultured bacterium Proteobacteria
OTU.212 11.743862 3.105905 0.6026082 4.739241 0.0000011 0.0000580 root_tissue Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingobium NA Proteobacteria
OTU.632 7.494015 3.038092 0.9908077 2.813959 0.0024468 0.0510916 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter Ambiguous_taxa Bacteroidetes
OTU.469 10.759062 3.028399 0.7534505 3.687566 0.0001132 0.0036174 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Niastella NA Bacteroidetes
OTU.1231 9.041788 3.005081 0.9467262 2.910114 0.0018065 0.0392770 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Flavobacteriaceae Flavobacterium Ambiguous_taxa Bacteroidetes
OTU.131 53.333470 2.969387 0.5904476 4.605637 0.0000021 0.0000923 root_tissue Bacteria Bacteroidetes Flavobacteriia Flavobacteriales Flavobacteriaceae Flavobacterium NA Bacteroidetes
OTU.1390 7.144050 2.959469 0.6741350 4.019180 0.0000292 0.0009775 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Ambiguous_taxa Ambiguous_taxa Proteobacteria
OTU.13254 4.714827 2.946382 0.9870268 2.731823 0.0031492 0.0638631 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter Ambiguous_taxa Bacteroidetes
OTU.1622 4.074434 2.940885 0.8745062 3.077034 0.0010454 0.0259372 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Niastella NA Bacteroidetes
OTU.93 5.181822 2.934541 0.9087872 2.953981 0.0015685 0.0974214 flowers Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Aurantimonadaceae Aureimonas Ambiguous_taxa Proteobacteria
OTU.8 67.362341 2.915148 0.3338322 7.983495 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium Proteobacteria
OTU.129 85.411783 2.911249 0.5686669 4.679803 0.0000014 0.0000704 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa Proteobacteria
OTU.15770 1.948971 2.879227 0.6513249 4.036736 0.0000271 0.0009368 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae NA NA Proteobacteria
OTU.16390 2.849999 2.857148 0.8100691 3.218427 0.0006445 0.0167803 root_tissue Bacteria Actinobacteria Actinobacteria Micromonosporales Micromonosporaceae NA NA Actinobacteria
OTU.552 11.016357 2.856283 0.6387420 4.080338 0.0000225 0.0008037 root_tissue Bacteria Proteobacteria Gammaproteobacteria Cellvibrionales Cellvibrionaceae Cellvibrio uncultured bacterium Proteobacteria
OTU.48 4.065260 2.855083 0.8595808 3.030644 0.0012202 0.0793429 flowers Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter Ambiguous_taxa Bacteroidetes
OTU.1348 5.812813 2.846243 1.0103038 2.569765 0.0050884 0.0949681 root_tissue Bacteria Proteobacteria Gammaproteobacteria Cellvibrionales Cellvibrionaceae Cellvibrio NA Proteobacteria
OTU.1186 5.167119 2.822710 0.7321251 3.514031 0.0002207 0.0066488 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae uncultured Ambiguous_taxa Bacteroidetes
OTU.7616 3.973050 2.817578 0.5908424 4.345622 0.0000069 0.0002712 root_tissue Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Aeromicrobium NA Actinobacteria
OTU.90 24.481151 2.770555 0.7447716 3.384333 0.0003568 0.0194346 leaves Bacteria Firmicutes Bacilli Bacillales Family XII Exiguobacterium NA Firmicutes
OTU.10124 4.065957 2.760471 0.6165990 4.071481 0.0000234 0.0008210 root_tissue Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces Ambiguous_taxa Actinobacteria
OTU.117 285.476713 2.752864 0.3483059 7.185821 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Aquabacterium NA Proteobacteria
OTU.7615 6.481045 2.718543 0.9334875 2.644431 0.0040914 0.0791634 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa Proteobacteria
OTU.74 149.983717 2.687463 0.4616853 5.279491 0.0000001 0.0000049 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium NA Proteobacteria
OTU.14360 2.249675 2.666589 0.7417945 3.257761 0.0005615 0.0151814 root_tissue Bacteria Proteobacteria Deltaproteobacteria Bdellovibrionales Bdellovibrionaceae Bdellovibrio NA Proteobacteria
OTU.8094 8.284363 2.631986 0.5230181 4.554309 0.0000026 0.0001131 root_tissue Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae NA NA Proteobacteria
OTU.224 6.274620 2.626845 0.8249601 2.881164 0.0019810 0.0953520 leaves Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa Proteobacteria
OTU.91 127.698326 2.622634 0.5217532 4.547425 0.0000027 0.0001145 root_tissue Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA Proteobacteria
OTU.31 142.057431 2.614664 0.8858850 2.669267 0.0038008 0.0742221 root_tissue Bacteria Acidobacteria Holophagae Holophagales Holophagaceae Holophaga NA Acidobacteria
OTU.548 12.116179 2.591488 0.8021953 2.918851 0.0017566 0.0385908 root_tissue Bacteria Proteobacteria Deltaproteobacteria Myxococcales BIrii41 uncultured bacterium uncultured bacterium Proteobacteria
OTU.16 132.514049 2.587894 0.4935809 4.736598 0.0000011 0.0001279 flowers Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA Actinobacteria
OTU.466 12.531213 2.570688 0.7589589 3.057726 0.0011151 0.0270320 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Chitinophaga NA Bacteroidetes
OTU.2470 4.613000 2.545342 0.7709293 2.977370 0.0014537 0.0329653 root_tissue Bacteria Verrucomicrobia Opitutae Opitutales Opitutaceae Opitutus uncultured bacterium Verrucomicrobia
OTU.388 39.532427 2.525511 0.7991196 2.847522 0.0022031 0.0469317 root_tissue Bacteria Chloroflexi Chloroflexia Herpetosiphonales Herpetosiphonaceae Herpetosiphon uncultured bacterium Chloroflexi
OTU.468 8.971922 2.516268 0.8488002 2.669966 0.0037929 0.0742221 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Sphingobacteriaceae Pedobacter NA Bacteroidetes
OTU.302 8.209868 2.493915 0.6576216 3.412168 0.0003222 0.0091839 root_tissue Bacteria Proteobacteria Deltaproteobacteria Oligoflexales 0319-6G20 uncultured bacterium uncultured bacterium Proteobacteria
OTU.4 28.819523 2.470702 0.6492221 3.420559 0.0003125 0.0177132 leaves Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Lactococcus Lactococcus lactis Firmicutes
OTU.819 16.916152 2.466713 0.5864983 3.779573 0.0000785 0.0025884 root_tissue Bacteria Verrucomicrobia Opitutae Opitutales Opitutaceae Opitutus NA Verrucomicrobia
OTU.49 477.208849 2.455047 0.3191560 6.908995 0.0000000 0.0000000 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae NA NA Proteobacteria
OTU.426 17.507103 2.422783 0.5036325 4.314223 0.0000080 0.0003071 root_tissue Bacteria Proteobacteria Betaproteobacteria Methylophilales Methylophilaceae Methylotenera NA Proteobacteria
OTU.280 16.454123 2.407027 0.7051562 3.058935 0.0011106 0.0270320 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae NA NA Proteobacteria
OTU.7137 6.963771 2.403990 0.5718004 3.767031 0.0000826 0.0026801 root_tissue Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Aeromicrobium uncultured bacterium Actinobacteria
OTU.9403 3.450851 2.378867 0.6601245 3.224947 0.0006300 0.0166078 root_tissue Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Nocardioides NA Actinobacteria
OTU.37 346.944734 2.358337 0.4175069 5.049826 0.0000002 0.0000161 leaves Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA Actinobacteria
OTU.624 4.904727 2.310261 0.6514538 3.162559 0.0007819 0.0198689 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA Proteobacteria
OTU.4825 9.177624 2.273039 0.6370144 3.175813 0.0007471 0.0192147 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Devosia NA Proteobacteria
OTU.11270 4.430509 2.210824 0.5979468 3.279262 0.0005204 0.0142534 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bosea uncultured bacterium Proteobacteria
OTU.9343 2.228053 2.201662 0.6665769 2.927887 0.0017064 0.0378814 root_tissue Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA Actinobacteria
OTU.59 103.911456 2.169586 0.7487657 2.563668 0.0051786 0.0958048 root_tissue Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Dyella Ambiguous_taxa Proteobacteria
OTU.120 6.328476 2.167536 0.6133107 3.126533 0.0008844 0.0607045 flowers Bacteria Actinobacteria Actinobacteria Kineosporiales Kineosporiaceae Kineococcus Ambiguous_taxa Actinobacteria
OTU.7 130.271901 2.161665 0.5720012 3.342064 0.0004158 0.0220975 leaves Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa Proteobacteria
OTU.36 124.294337 2.141622 0.3952198 4.786253 0.0000008 0.0000472 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Devosia NA Proteobacteria
OTU.145 55.194300 2.105192 0.4134932 4.486631 0.0000036 0.0001469 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Phyllobacteriaceae Mesorhizobium NA Proteobacteria
OTU.1005 3.305019 2.092650 0.6215192 2.964751 0.0015146 0.0339827 root_tissue Bacteria Actinobacteria Acidimicrobiia Acidimicrobiales uncultured NA NA Actinobacteria
OTU.6 1346.830825 2.055248 0.5600897 3.223141 0.0006340 0.0474706 flowers Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Bacillus cereus Firmicutes
OTU.13 238.573815 2.016071 0.3238031 5.454152 0.0000000 0.0000032 flowers Bacteria Proteobacteria Betaproteobacteria Burkholderiales Burkholderiaceae Ralstonia NA Proteobacteria
OTU.366 19.928246 1.981334 0.4721376 3.667010 0.0001227 0.0038623 root_tissue Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Nocardioides uncultured bacterium Actinobacteria
OTU.571 22.623607 1.872554 0.5366803 3.023315 0.0012501 0.0296233 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa NA Proteobacteria
OTU.672 15.946577 1.863494 0.4446025 3.629071 0.0001422 0.0043470 root_tissue Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Novosphingobium NA Proteobacteria
OTU.2104 9.875240 1.861001 0.6313339 2.551743 0.0053593 0.0975347 root_tissue Bacteria Bacteroidetes Sphingobacteriia Sphingobacteriales Chitinophagaceae Flavisolibacter uncultured bacterium Bacteroidetes
OTU.2341 3.022184 1.848761 0.5641277 2.834040 0.0022982 0.0484685 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae NA NA Proteobacteria
OTU.1116 2.617411 1.827161 0.6066193 2.599919 0.0046623 0.0877925 root_tissue Bacteria Firmicutes Bacilli Bacillales Paenibacillaceae Paenibacillus NA Firmicutes
OTU.743 12.443579 1.793824 0.6050902 2.551394 0.0053646 0.0975347 root_tissue Bacteria Proteobacteria Deltaproteobacteria Myxococcales Haliangiaceae Haliangium Ambiguous_taxa Proteobacteria
OTU.9296 104.815616 1.634990 0.4595231 3.013972 0.0012893 0.0302116 root_tissue Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA Actinobacteria
OTU.310 16.316297 1.575165 0.5090802 2.603058 0.0046198 0.0877925 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Phyllobacteriaceae Mesorhizobium Mesorhizobium ciceri Proteobacteria
OTU.50 18.537429 1.573867 0.4924443 2.688359 0.0035902 0.0714316 root_tissue Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium NA Proteobacteria
OTU.384 79.199405 1.530956 0.3883509 3.298449 0.0004861 0.0134893 root_tissue Bacteria Actinobacteria Actinobacteria Propionibacteriales Nocardioidaceae Aeromicrobium Ambiguous_taxa Actinobacteria
OTU.18 313.470752 1.525645 0.2856414 4.465897 0.0000040 0.0001586 root_tissue Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA Proteobacteria
OTU.61 93.231576 1.519089 0.4074767 3.114508 0.0009213 0.0615253 flowers Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Ambiguous_taxa Firmicutes
OTU.79 8.780912 1.489666 0.2869787 4.319715 0.0000078 0.0605239 rhizosphere_soil Bacteria Actinobacteria Actinobacteria Pseudonocardiales Pseudonocardiaceae Ambiguous_taxa Ambiguous_taxa Actinobacteria
OTU.19 368.944068 1.464972 0.3683052 3.298819 0.0004855 0.0134893 root_tissue Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA Actinobacteria
OTU.140 55.753479 1.274913 0.3940312 2.601095 0.0046463 0.0877925 root_tissue Bacteria Proteobacteria Alphaproteobacteria Caulobacterales Caulobacteraceae Phenylobacterium uncultured bacterium Proteobacteria

Fungi

fungi.deseq.df = data.frame()
for (pcomp in c("rhizosphere_soil", "leaves", "flowers")){
  # Subset the data to just include bulk soil and plant compartment samples.
  physeqITS.fungi.BSPC = subset_samples(physeqITS.fungi, Sample_type %in% c("bulk_soil", pcomp))
  physeqITS.fungi.BSPC = prune_taxa(taxa_sums(physeqITS.fungi.BSPC) > 0, physeqITS.fungi.BSPC)
  
  fungi.metadata = data.frame(sample_data(physeqITS.fungi.BSPC), stringsAsFactors = F)
  
  # Set the sample types as a factor
  sample_data(physeqITS.fungi.BSPC)$Sample_type = factor(sample_data(physeqITS.fungi.BSPC)$Sample_type, levels = c("bulk_soil", pcomp))
  sample_data(physeqITS.fungi.BSPC)$Field = as.factor(sample_data(physeqITS.fungi.BSPC)$Field)
  
  # Get the taxonomy table
  fungi.tax.table = data.frame(gsub("D_.__", "", tax_table(physeqITS.fungi.BSPC)), stringsAsFactors = F) %>%
    rownames_to_column(var="OTU") %>%
    dplyr::select(OTU, Kingdom, Phylum, Class, Order, Family, Genus, Species)
  
  # I'll also remove sparse taxa. In this case, taxa must be found with at least a count of 5 in at least 10 samples.
  OTU.table = as.matrix(otu_table(physeqITS.fungi.BSPC))
  OTU.table[OTU.table < 5] = 0
  OTU.table[OTU.table >= 5] = 1
  #sparse = length(sample_names(physeqITS.fungi.BSPC)) * 0.25
  sparse = 1
  sampleCount.df = data.frame(sample = rowSums(OTU.table)) %>%
    rownames_to_column(var="OTU") %>%
    filter(sample >= sparse)
  physeqITS.fungi.BSPC.sparse = prune_taxa(sampleCount.df$OTU, physeqITS.fungi.BSPC)
  
  # Run DESeq2
  fungi.deseq = phyloseq_to_deseq2(physeqITS.fungi.BSPC.sparse, ~ Field + Sample_type)
  fungi.deseq = DESeq(fungi.deseq, betaPrior=TRUE)
  fungi.deseq.res = results(fungi.deseq, lfcThreshold = .25,
                      contrast = c("Sample_type", pcomp, "bulk_soil"), 
                      altHypothesis = "greater",
                      test="Wald")
  fungi.deseq.df = rbind(fungi.deseq.df, data.frame(fungi.deseq.res) %>%
                          rownames_to_column(var="OTU") %>%
                          mutate(Sample_type = pcomp))
}

fungi.deseq.df = inner_join(fungi.deseq.df, fungi.tax.table, by="OTU") %>%
  filter(!(is.na(padj))) %>%
  arrange(-log2FoldChange) %>%
  mutate(sig = ifelse(padj <= 0.1, Phylum, "not significant"))
fungi.deseq.df$Phylum = factor(fungi.deseq.df$Phylum, levels = unique(fungi.deseq.df$Phylum))

kable(fungi.deseq.df %>% filter(padj < 0.1)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
OTU baseMean log2FoldChange lfcSE stat pvalue padj Sample_type Kingdom Phylum Class Order Family Genus Species sig
OTU_3520 312.894261 6.560139 2.3833162 2.647630 0.0040529 0.0333352 leaves Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU Basidiomycota
OTU_27 738.541513 6.418929 0.4133698 14.923510 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU Basidiomycota
OTU_331 41.448986 6.377231 0.5055712 12.119423 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Dioszegia Dioszegia_hungarica_SH182099.07FU Basidiomycota
OTU_338 34.257294 6.213597 0.5261222 11.335003 0.0000000 0.0000000 leaves Fungi Basidiomycota NA NA NA NA NA Basidiomycota
OTU_336 36.932173 5.918484 0.6866092 8.255764 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Vishniacozyma Vishniacozyma_victoriae_SH181628.07FU Basidiomycota
OTU_118 67.337604 5.905907 0.4972883 11.373495 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales Pleosporaceae Alternaria Alternaria_infectoria_SH434793.07FU Ascomycota
OTU_21 1960.735503 5.879615 0.4239653 13.278480 0.0000000 0.0000000 leaves Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU Basidiomycota
OTU_405 19.166518 5.764982 0.5112402 10.787456 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales NA NA NA Ascomycota
OTU_84 98.385896 5.718279 0.6134294 8.914276 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Neoascochyta NA Ascomycota
OTU_333 44.967451 5.680707 0.5134989 10.575889 0.0000000 0.0000000 leaves Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU Basidiomycota
OTU_346 24.852865 5.584749 0.4735202 11.266148 0.0000000 0.0000000 leaves Fungi NA NA NA NA NA NA NA
OTU_174 812.817534 5.580714 0.4152821 12.836370 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA Ascomycota
OTU_359 34.849124 5.448842 0.5438816 9.558777 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales NA NA NA Ascomycota
OTU_373 16.542217 5.282998 0.6397823 7.866736 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_crocea_SH209709.07FU Basidiomycota
OTU_360 17.538441 5.013087 0.5169388 9.214025 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Capnodiales Mycosphaerellaceae NA NA Ascomycota
OTU_343 26.461883 4.976994 0.4969621 9.511780 0.0000000 0.0000000 leaves Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_roseus_SH194973.07FU Basidiomycota
OTU_259 14.978191 4.835848 0.7761212 5.908675 0.0000000 0.0000000 leaves Fungi NA NA NA NA NA NA NA
OTU_375 13.466806 4.364706 0.6468391 6.361251 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Bulleribasidium Bulleribasidium_oberjochense_SH222253.07FU Basidiomycota
OTU_1247 10.001585 4.207004 0.5996886 6.598432 0.0000000 0.0000000 leaves Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU Basidiomycota
OTU_251 10.652648 3.934823 0.6300000 5.848926 0.0000000 0.0000000 leaves Fungi Ascomycota Dothideomycetes Pleosporales NA NA NA Ascomycota
OTU_232 11.359962 3.795625 0.5363343 6.610848 0.0000000 0.0000000 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Hannaella Hannaella_coprosmae_SH219462.07FU Basidiomycota
OTU_402 7.773317 3.609888 0.7768567 4.324978 0.0000076 0.0000965 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Dioszegia Dioszegia_changbaiensis_SH209715.07FU Basidiomycota
OTU_340 13.029190 3.604033 0.6983277 4.802949 0.0000008 0.0000107 leaves Fungi Basidiomycota Tremellomycetes Tremellales Rhynchogastremataceae Papiliotrema Papiliotrema_fuscus_SH205049.07FU Basidiomycota
OTU_332 16.394310 3.557134 0.6300355 5.249123 0.0000001 0.0000011 leaves Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium NA Basidiomycota
OTU_350 7.029060 3.386908 0.6756842 4.642566 0.0000017 0.0000226 leaves Fungi Ascomycota Taphrinomycetes Taphrinales Taphrinaceae Taphrina NA Ascomycota
OTU_125 187.638569 3.371978 0.4444258 7.024746 0.0000000 0.0000000 leaves Fungi NA NA NA NA NA NA NA
OTU_736 6.328401 3.125325 0.7598375 3.784131 0.0000771 0.0008458 leaves Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium Filobasidium_floriforme_SH264587.07FU Basidiomycota
OTU_337 9.583233 3.123467 0.6844607 4.198148 0.0000135 0.0001640 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Vishniacozyma Vishniacozyma_carnescens_SH181629.07FU Basidiomycota
OTU_616 3.265198 2.926202 0.8783275 3.046929 0.0011560 0.0111856 leaves Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium NA Basidiomycota
OTU_339 9.533473 2.858591 0.9213092 2.831396 0.0023173 0.0206049 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Vishniacozyma Vishniacozyma_victoriae_SH181632.07FU Basidiomycota
OTU_21 143.719806 2.853298 0.5018934 5.186954 0.0000001 0.0000265 rhizosphere_soil Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU Basidiomycota
OTU_404 16.196603 2.799823 0.6399170 3.984615 0.0000338 0.0003834 leaves Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Neoascochyta Neoascochyta_desmazieri_SH215819.07FU Ascomycota
OTU_376 7.560157 2.784112 0.6920096 3.661961 0.0001251 0.0013282 leaves Fungi Basidiomycota Spiculogloeomycetes Spiculogloeales Spiculogloeaceae Phyllozyma NA Basidiomycota
OTU_333 6.711295 2.533756 0.6790831 3.363000 0.0003855 0.0929060 flowers Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU Basidiomycota
OTU_428 4.490908 2.498277 0.8609669 2.611340 0.0045094 0.0361852 leaves Fungi Basidiomycota Tremellomycetes Tremellales NA NA NA Basidiomycota
OTU_1902 2.736154 2.476924 0.7587722 2.934904 0.0016683 0.0152460 leaves Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae NA NA Ascomycota
OTU_605 5.928250 2.450875 0.5474040 4.020567 0.0000290 0.0003411 leaves Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_patagonicus_SH194975.07FU Basidiomycota
OTU_27 16.824086 2.432515 0.4146040 5.264097 0.0000001 0.0000265 rhizosphere_soil Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU Basidiomycota
OTU_363 4.449375 2.352861 0.6578100 3.196760 0.0006949 0.0069280 leaves Fungi Basidiomycota Tremellomycetes Tremellales Tremellaceae Cryptococcus Cryptococcus_perniciosus_SH205048.07FU Basidiomycota
OTU_430 2.554322 2.307994 0.8427208 2.442083 0.0073014 0.0533814 leaves Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium NA Basidiomycota
OTU_425 2.705111 2.290928 0.6063100 3.366146 0.0003811 0.0039185 leaves Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_symmetricus_SH205939.07FU Basidiomycota
OTU_408 4.391839 2.235298 0.6762992 2.935533 0.0016649 0.0152460 leaves Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium Filobasidium_oeirense_SH264590.07FU Basidiomycota
OTU_604 3.382446 2.231026 0.7350034 2.695261 0.0035167 0.0296663 leaves Fungi NA NA NA NA NA NA NA
OTU_115 252.395579 2.089184 0.2542654 7.233324 0.0000000 0.0000000 leaves Fungi Ascomycota NA NA NA NA NA Ascomycota
OTU_27 87.776826 2.067084 0.4878564 3.724629 0.0000978 0.0471402 flowers Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU Basidiomycota
OTU_611 2.078456 2.038240 0.7203751 2.482374 0.0065255 0.0487931 leaves Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Vishniacozyma Vishniacozyma_dimennae_SH203798.07FU Basidiomycota
OTU_594 1.507002 2.023885 0.6519145 2.721039 0.0032539 0.0281715 leaves Fungi Ascomycota Dothideomycetes Pleosporales Massarinaceae Stagonospora Stagonospora_pseudovitensis_SH182977.07FU Ascomycota
OTU_386 2.415802 1.958120 0.7794132 2.191546 0.0142061 0.0994430 leaves Fungi NA NA NA NA NA NA NA
OTU_345 64.780904 1.953952 0.6714589 2.537687 0.0055794 0.0426888 leaves Fungi Basidiomycota Tremellomycetes Tremellales Trimorphomycetaceae Saitozyma Saitozyma_paraflava_SH181618.07FU Basidiomycota
OTU_378 2.770201 1.882161 0.6889872 2.368928 0.0089199 0.0637964 leaves Fungi Ascomycota Dothideomycetes Capnodiales Dissoconiaceae Dissoconium NA Ascomycota
OTU_422 18.178322 1.616271 0.5311331 2.572370 0.0050502 0.0395602 leaves Fungi Ascomycota Dothideomycetes Pleosporales Didymosphaeriaceae Pseudopithomyces Pseudopithomyces_chartarum_SH630197.07FU Ascomycota

Part 4: Core OTUs

For defining the core microbiome we can find OTUs that are present on all plants. Again this might mean that they are just super abundant in the environment though. Then we can narrow this list by seeing which OTUs are enriched on the plant compartment compared to the bulk soil. This would suggest that, while they may be found in the soil, as this might be the source of the microbiome, they may be more associated with the plant.

In part 2 and 3 we found the ubiquitous OTUs and the significantly enriched OTUs respectively. So lets now merge the data together to identify the Core OTUs. Also lets still include the top 5 most abundant OTUs.

Bacteria

## Get the ubiqitous OTUs and whether they are CORE or not (significantly enriched)
bact.ubiq = bact.OTU.ubiq %>%
  mutate(Ubiq_in_plant = "YES") %>%
  dplyr::select(Sample_type, OTU, Ubiq_in_plant)
bact.enr = bact.deseq.df %>% 
  dplyr::select(Sample_type, OTU, log2FoldChange, lfcSE, padj)
  
## Whether OTUS are ubiquitous in the bulk soil
bact.OTU.soil.ubiq = bact.OTU.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  mutate(Ubiq_in_soil = ifelse(min_RA > 0, "YES", "NO")) %>%
  dplyr::select(OTU, Ubiq_in_soil)

## Final table of OTUS of interest including CORE
bact.OTU.interest = bact.OTU.rare.sum %>%
  arrange(-mean_RA) %>%
  group_by(Sample_type) %>%
  mutate(RA_Rank = row_number()) %>%
  ungroup() %>%
  dplyr::select(Sample_type, OTU, mean_RA, SE_RA, RA_Rank, Rank1, Rank2, Rank3, Rank4, Rank5, Rank6, Rank7) %>%
  full_join(bact.OTU.soil.ubiq, by="OTU") %>%
  full_join(bact.ubiq, by=c("Sample_type", "OTU")) %>%
  full_join(bact.enr, by=c("Sample_type", "OTU")) %>%
  mutate(Ubiq_in_plant = ifelse(is.na(Ubiq_in_plant), "NO", Ubiq_in_plant),
         Core = ifelse(Ubiq_in_plant == "YES" & padj < 0.1, "YES", "NO"),
         mean_RA = mean_RA*100,
         SE_RA = SE_RA*100,
         log2FoldChange = ifelse(padj > 0.1, NA, log2FoldChange),
         lfcSE = ifelse(padj > 0.1, NA, lfcSE)) %>%
  filter(RA_Rank <= 5 | Core == "YES") %>%
  dplyr::select(Sample_type, OTU, mean_RA, SE_RA, log2FoldChange, lfcSE, RA_Rank, Ubiq_in_plant, Ubiq_in_soil, Core, Rank1, Rank2, Rank3, Rank4, Rank5, Rank6, Rank7) %>%
  dplyr::rename(Domain=Rank1, Phylum=Rank2, Class=Rank3, Order=Rank4, Family=Rank5, Genus=Rank6, Species=Rank7) %>%
  mutate(Sample_type = factor(Sample_type, levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"))) %>%
  arrange(Sample_type, RA_Rank)

kable(bact.OTU.interest %>%
        mutate(RA = paste(round(mean_RA, digits = 2), " (", round(SE_RA, digits = 2), ")", sep=""),
        log2FC = paste(round(log2FoldChange, digits = 2), " (", round(lfcSE, digits = 2), ")", sep="")) %>%
  dplyr::select(Sample_type, OTU, RA, log2FC, RA_Rank, Ubiq_in_plant, Ubiq_in_soil, Core, Domain, Phylum, Class, Order, Family, Genus, Species)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
Sample_type OTU RA log2FC RA_Rank Ubiq_in_plant Ubiq_in_soil Core Domain Phylum Class Order Family Genus Species
bulk_soil OTU.10 1.77 (0.25) NA (NA) 1 NO YES NO Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Pseudarthrobacter NA
bulk_soil OTU.28 1.67 (0.14) NA (NA) 2 NO NO NO Bacteria Chloroflexi KD4-96 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
bulk_soil OTU.21 1.29 (0.13) NA (NA) 3 NO YES NO Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Janibacter NA
bulk_soil OTU.66 1.23 (0.23) NA (NA) 4 NO YES NO Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
bulk_soil OTU.39 1.12 (0.14) NA (NA) 5 NO YES NO Bacteria Proteobacteria Gammaproteobacteria Xanthomonadales Xanthomonadaceae Lysobacter NA
rhizosphere_soil OTU.10 2 (0.22) NA (NA) 1 YES YES NO Bacteria Actinobacteria Actinobacteria Micrococcales Micrococcaceae Pseudarthrobacter NA
rhizosphere_soil OTU.28 1.63 (0.12) NA (NA) 2 NO NO NO Bacteria Chloroflexi KD4-96 uncultured bacterium uncultured bacterium uncultured bacterium uncultured bacterium
rhizosphere_soil OTU.21 1.2 (0.16) NA (NA) 3 NO YES NO Bacteria Actinobacteria Actinobacteria Micrococcales Intrasporangiaceae Janibacter NA
rhizosphere_soil OTU.66 1 (0.11) NA (NA) 4 YES YES NO Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Ambiguous_taxa Ambiguous_taxa
rhizosphere_soil OTU.63 0.85 (0.07) NA (NA) 5 NO NO NO Bacteria Acidobacteria Subgroup 6 NA NA NA NA
root_tissue OTU.24 5.69 (0.79) 4.1 (0.38) 1 NO NO NO Bacteria Actinobacteria Actinobacteria Streptomycetales Streptomycetaceae Streptomyces Ambiguous_taxa
root_tissue OTU.17 4.94 (1.33) 8.62 (0.72) 2 NO NO NO Bacteria Proteobacteria Gammaproteobacteria Cellvibrionales Halieaceae NA NA
root_tissue OTU.49 3.47 (0.65) 2.46 (0.32) 3 YES YES YES Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae NA NA
root_tissue OTU.18 2.68 (0.6) 1.53 (0.29) 4 YES YES YES Bacteria Proteobacteria Betaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
root_tissue OTU.117 2.43 (0.45) 2.75 (0.35) 5 YES YES YES Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Aquabacterium NA
root_tissue OTU.74 1.09 (0.23) 2.69 (0.46) 14 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Rhizobiaceae Rhizobium NA
root_tissue OTU.91 1.05 (0.35) 2.62 (0.52) 17 YES NO YES Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
root_tissue OTU.8 0.59 (0.12) 2.92 (0.33) 32 YES YES YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
leaves OTU.8 13.14 (1.88) 6.69 (0.39) 1 YES YES YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
leaves OTU.5 12.81 (1.99) 7.65 (0.46) 2 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
leaves OTU.9 8.55 (1.73) 7.46 (0.57) 3 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa
leaves OTU.11 6.53 (1.09) 7.17 (0.59) 4 YES NO YES Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves OTU.23 5.63 (1.78) 8.3 (0.67) 5 YES NO YES Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves OTU.16 4.91 (0.83) 4.98 (0.44) 6 YES YES YES Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
leaves OTU.37 2.28 (0.36) 2.36 (0.42) 8 YES YES YES Bacteria Actinobacteria Actinobacteria Micrococcales Microbacteriaceae NA NA
leaves OTU.25 1.67 (0.45) 7.53 (0.75) 9 YES NO YES Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Hymenobacter uncultured bacterium
leaves OTU.3 1.53 (0.34) 5.57 (0.61) 11 YES NO YES Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
leaves OTU.137 0.49 (0.09) 6.47 (0.58) 30 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
leaves OTU.1896 0.23 (0.04) 4.6 (0.53) 44 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA
flowers OTU.7 12.56 (3.19) 6.21 (0.59) 1 NO NO NO Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas Ambiguous_taxa
flowers OTU.4 11.67 (4.41) 9.92 (0.72) 2 NO NO NO Bacteria Firmicutes Bacilli Lactobacillales Streptococcaceae Lactococcus Lactococcus lactis
flowers OTU.12 6.33 (2.01) 8.03 (0.73) 3 NO NO NO Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Pantoea Ambiguous_taxa
flowers OTU.20 6.23 (2) 5.85 (0.62) 4 NO NO NO Bacteria Proteobacteria Gammaproteobacteria Enterobacteriales Enterobacteriaceae Enterobacter NA
flowers OTU.6 5.53 (3.59) 2.06 (0.56) 5 YES YES YES Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Bacillus cereus
flowers OTU.13 2.85 (0.85) 2.02 (0.32) 7 YES YES YES Bacteria Proteobacteria Betaproteobacteria Burkholderiales Burkholderiaceae Ralstonia NA
flowers OTU.61 1.22 (0.55) 1.52 (0.41) 10 YES YES YES Bacteria Firmicutes Bacilli Bacillales Bacillaceae Bacillus Ambiguous_taxa
flowers OTU.5 1.05 (0.36) 4.94 (0.55) 11 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas Ambiguous_taxa
flowers OTU.8 1.01 (0.33) 4.07 (0.46) 12 YES YES YES Bacteria Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas uncultured bacterium
flowers OTU.9568 0.72 (0.18) 6.74 (0.56) 16 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Bradyrhizobiaceae Bradyrhizobium NA
flowers OTU.9 0.37 (0.13) 3.41 (0.54) 24 YES NO YES Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Methylobacteriaceae Methylobacterium Ambiguous_taxa

Fungi

## Get the ubiqitous OTUs and whether they are CORE or not (significantly enriched)
fungi.ubiq = fungi.OTU.ubiq %>%
  mutate(Ubiq_in_plant = "YES") %>%
  dplyr::select(Sample_type, OTU, Ubiq_in_plant)
fungi.enr = fungi.deseq.df %>% 
  dplyr::select(Sample_type, OTU, log2FoldChange, lfcSE, padj)
  
## Whether OTUS are ubiquitous in the bulk soil
fungi.OTU.soil.ubiq = fungi.OTU.sum %>%
  filter(Sample_type == "bulk_soil") %>%
  mutate(Ubiq_in_soil = ifelse(min_RA > 0, "YES", "NO")) %>%
  dplyr::select(OTU, Ubiq_in_soil)

## Final table of OTUS of interest including CORE
fungi.OTU.interest = fungi.OTU.rare.sum %>%
  arrange(-mean_RA) %>%
  group_by(Sample_type) %>%
  mutate(RA_Rank = row_number()) %>%
  ungroup() %>%
  dplyr::select(Sample_type, OTU, mean_RA, SE_RA, RA_Rank, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
  full_join(fungi.OTU.soil.ubiq, by="OTU") %>%
  full_join(fungi.ubiq, by=c("Sample_type", "OTU")) %>%
  full_join(fungi.enr, by=c("Sample_type", "OTU")) %>%
  filter(RA_Rank <= 5 | Ubiq_in_plant == "YES") %>%
   mutate(Ubiq_in_plant = ifelse(is.na(Ubiq_in_plant), "NO", Ubiq_in_plant),
         Core = ifelse(Ubiq_in_plant == "YES" & padj < 0.1, "YES", "NO"),
         mean_RA = mean_RA*100,
         SE_RA = SE_RA*100,
         log2FoldChange = ifelse(padj > 0.1, NA, log2FoldChange),
         lfcSE = ifelse(padj > 0.1, NA, lfcSE)) %>%
  filter(RA_Rank <= 5 | Core == "YES") %>%
  dplyr::select(Sample_type, OTU, mean_RA, SE_RA, log2FoldChange, lfcSE, RA_Rank, Ubiq_in_plant, Ubiq_in_soil, Core, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
  mutate(Sample_type = factor(Sample_type, levels = c("bulk_soil", "rhizosphere_soil", "root_tissue", "leaves", "flowers"))) %>%
  arrange(Sample_type, RA_Rank)

kable(fungi.OTU.interest %>%
        mutate(RA = paste(round(mean_RA, digits = 2), " (", round(SE_RA, digits = 2), ")", sep=""),
               log2FC = paste(round(log2FoldChange, digits = 2), " (", round(lfcSE, digits = 2), ")", sep="")) %>%
  dplyr::select(Sample_type, OTU, RA, log2FC, RA_Rank, Ubiq_in_plant, Ubiq_in_soil, Core, Kingdom, Phylum, Class, Order, Family, Genus, Species)) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height="400px")
Sample_type OTU RA log2FC RA_Rank Ubiq_in_plant Ubiq_in_soil Core Kingdom Phylum Class Order Family Genus Species
bulk_soil OTU_13 13.97 (2.28) NA (NA) 1 NO YES NO Fungi Ascomycota Sordariomycetes Glomerellales Plectosphaerellaceae Verticillium Verticillium_dahliae_SH466940.07FU
bulk_soil OTU_8 6.93 (0.98) NA (NA) 2 NO YES NO Fungi Ascomycota Sordariomycetes Hypocreales Nectriaceae Fusarium NA
bulk_soil OTU_34 6.84 (1.5) NA (NA) 3 NO YES NO Fungi NA NA NA NA NA NA
bulk_soil OTU_309 5.08 (1.52) NA (NA) 4 NO YES NO Fungi NA NA NA NA NA NA
bulk_soil OTU_33 4.79 (2.07) NA (NA) 5 NO YES NO Fungi Ascomycota Sordariomycetes Hypocreales Hypocreales_fam_Incertae_sedis NA NA
rhizosphere_soil OTU_13 10.94 (1.6) NA (NA) 1 YES YES NO Fungi Ascomycota Sordariomycetes Glomerellales Plectosphaerellaceae Verticillium Verticillium_dahliae_SH466940.07FU
rhizosphere_soil OTU_34 9.57 (1.81) NA (NA) 2 YES YES NO Fungi NA NA NA NA NA NA
rhizosphere_soil OTU_8 7.63 (1.57) NA (NA) 3 YES YES NO Fungi Ascomycota Sordariomycetes Hypocreales Nectriaceae Fusarium NA
rhizosphere_soil OTU_21 4.66 (1.45) 2.85 (0.5) 4 NO NO NO Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
rhizosphere_soil OTU_38 4.31 (0.69) NA (NA) 5 YES YES NO Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Didymella Didymella_dimorpha_SH245093.07FU
rhizosphere_soil OTU_27 0.46 (0.18) 2.43 (0.41) 46 YES NO YES Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
leaves OTU_21 36.25 (3.11) 5.88 (0.42) 1 YES NO YES Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
leaves OTU_174 16.75 (1.67) 5.58 (0.42) 2 YES NO YES Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
leaves OTU_27 11.64 (1.64) 6.42 (0.41) 3 YES NO YES Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
leaves OTU_115 5.27 (0.68) 2.09 (0.25) 4 YES YES YES Fungi Ascomycota NA NA NA NA NA
leaves OTU_125 4.15 (0.59) 3.37 (0.44) 5 YES NO YES Fungi NA NA NA NA NA NA
leaves OTU_84 2.21 (0.44) 5.72 (0.61) 8 YES NO YES Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Neoascochyta NA
leaves OTU_118 1.74 (0.33) 5.91 (0.5) 9 YES NO YES Fungi Ascomycota Dothideomycetes Pleosporales Pleosporaceae Alternaria Alternaria_infectoria_SH434793.07FU
leaves OTU_331 0.88 (0.12) 6.38 (0.51) 10 YES NO YES Fungi Basidiomycota Tremellomycetes Tremellales Bulleribasidiaceae Dioszegia Dioszegia_hungarica_SH182099.07FU
leaves OTU_338 0.8 (0.12) 6.21 (0.53) 11 YES NO YES Fungi Basidiomycota NA NA NA NA NA
leaves OTU_333 0.75 (0.16) 5.68 (0.51) 12 YES NO YES Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU
leaves OTU_346 0.7 (0.13) 5.58 (0.47) 13 YES NO YES Fungi NA NA NA NA NA NA
leaves OTU_359 0.6 (0.07) 5.45 (0.54) 16 YES NO YES Fungi Ascomycota Dothideomycetes Pleosporales NA NA NA
leaves OTU_332 0.36 (0.08) 3.56 (0.63) 20 YES NO YES Fungi Basidiomycota Tremellomycetes Filobasidiales Filobasidiaceae Filobasidium NA
leaves OTU_360 0.27 (0.04) 5.01 (0.52) 27 YES NO YES Fungi Ascomycota Dothideomycetes Capnodiales Mycosphaerellaceae NA NA
flowers OTU_21 32.3 (5.16) NA (NA) 1 YES NO NO Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
flowers OTU_115 22.33 (3.35) NA (NA) 2 YES YES NO Fungi Ascomycota NA NA NA NA NA
flowers OTU_27 10.8 (2.24) 2.07 (0.49) 3 YES NO YES Fungi Basidiomycota Tremellomycetes Tremellales Bulleraceae Bullera Bullera_alba_SH215453.07FU
flowers OTU_174 8.18 (1.99) NA (NA) 4 YES NO NO Fungi Ascomycota Dothideomycetes Pleosporales Didymellaceae Epicoccum NA
flowers OTU_3520 3.22 (2.26) NA (NA) 5 NO NO NO Fungi Basidiomycota Exobasidiomycetes Entylomatales Entylomatales_fam_Incertae_sedis Tilletiopsis Tilletiopsis_washingtonensis_SH186666.07FU
flowers OTU_333 1.43 (0.31) 2.53 (0.68) 8 YES NO YES Fungi Basidiomycota Microbotryomycetes Sporidiobolales Sporidiobolaceae Sporobolomyces Sporobolomyces_ruberrimus_SH194976.07FU

Part 5: Crop yeild correlations

While we dont have yeild measures for each individual plants we do have them for the entire field. To see if there is any connection between the microbiomes and the hemp yeild lets try to find a correlation between the difference in grain yeild between fields and the community distance (Bray-Curtis).

Bacteria

for (plant_part in c("rhizosphere_soil", "root_tissue", "leaves", "flowers")){
  writeLines(paste("Running", plant_part, "\n", sep=" "))
  # Subset phyloseq
  physeq.sub = subset_samples(physeq16S.bact, Sample_type == plant_part & Field != "Freeville")
  
  # Remove fields with fewer than 3 replicates
  rep.n = data.frame(sample_data(physeq.sub), stringsAsFactors = F) %>%
    group_by(Field) %>%
    mutate(N = n()) %>%
    ungroup
  physeq.sub = prune_samples(rep.n[rep.n$N >= 3,]$sample_ID, physeq.sub)
  writeLines(paste("Removed fields:", unique(rep.n[rep.n$N < 3,]$Field)))
  
  # rarefy and normalize counts
  set.seed(72)
  physeq.sub = rarefy_even_depth(physeq.sub)
  physeq.sub = transform_sample_counts(physeq.sub, function(x) x / sum(x))
  
  # get sample metadata and split into 2 similar dataframes for contrasts
  metadata = data.frame(sample_data(physeq.sub), stringsAsFactors = F)
  meta1 = metadata %>%
    dplyr::select(sample_ID, Field)
  colnames(meta1) = c("sample1", "field1")
  meta2 = metadata %>%
    dplyr::select(sample_ID, Field)
  colnames(meta2) = c("sample2", "field2")
  
  # Calculate pairwise Bray-Curtis distances and add metadata
  BC.dist = data.frame(as.matrix(phyloseq::distance(physeq.sub, method="bray")))
  colnames(BC.dist) = gsub("X", "", gsub("\\.", "-", colnames(BC.dist)))
  dist.df = BC.dist %>%
    rownames_to_column(var="sample1") %>%
    gather(key="sample2", value="BC_dist", -sample1) %>%
    filter(sample1 != sample2) %>%
    mutate(contrast = paste(pmin(sample1, sample2), pmax(sample1, sample2), sep="+")) %>%
    separate(contrast, into=c("sample1", "sample2"), sep="\\+", remove=F) %>%
    unique() %>%
    inner_join(meta1) %>%
    inner_join(meta2) %>%
    filter(field1 != field2) %>%
    group_by(field1, field2) %>%
    summarize(mean_BC = mean(BC_dist),
              sd_BC = sd(BC_dist),
              n.con = n())
  
  # Add field data
  field.diff1 = field.data %>%
    dplyr::select(Field, Clean_Grain_Yield)
  colnames(field.diff1) = c("field1", "CGY1")
  field.diff2 = field.data %>%
    dplyr::select(Field, Clean_Grain_Yield)
  colnames(field.diff2) = c("field2", "CGY2")
  
  dist.field.df = inner_join(dist.df, field.diff1, by="field1")
  dist.field.df = inner_join(dist.field.df, field.diff2, by="field2")
  dist.field.df = dist.field.df %>%
    mutate(CGY_diff = abs(CGY1-CGY2))
  
  # Calculate correlation coeficient
  writeLines("Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield")
  corr = cor.test(dist.field.df$mean_BC, dist.field.df$CGY_diff, method = "pearson", conf.level = 0.95)
  print(corr)
  writeLines(paste("Adjusted pvalue =", p.adjust(corr$p.value, method="bonferroni", n=4), "\n"))
  
  # Plot
  title = paste(plant_part, ":", " r=", round(corr$estimate, digits=5), ", df=", corr$parameter, ", p=", round(p.adjust(corr$p.value, method="bonferroni", n=4), digits=5), sep="")
  dist.field.plot = ggplot(data=dist.field.df, aes(x=mean_BC, y=CGY_diff)) +
    geom_point() +
    geom_smooth(method="lm") +
    labs(x="Mean Bray-Curtis dissimilarity", y="Difference in clean grain yield") +
    ggtitle(label = title) +
    theme_bw()
  
  plot(dist.field.plot)
  
}
## Running rhizosphere_soil 
## 
## Removed fields: 
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = -0.30047, df = 8, p-value = 0.7715
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6894098  0.5613235
## sample estimates:
##        cor 
## -0.1056377 
## 
## Adjusted pvalue = 1

## Running root_tissue 
## 
## Removed fields: 
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = 0.76451, df = 8, p-value = 0.4665
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4411741  0.7648944
## sample estimates:
##       cor 
## 0.2609327 
## 
## Adjusted pvalue = 1

## Running leaves 
## 
## Removed fields: 
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = -1.9574, df = 8, p-value = 0.086
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.88249688  0.09436716
## sample estimates:
##        cor 
## -0.5690711 
## 
## Adjusted pvalue = 0.34399335864368

## Running flowers 
## 
## Removed fields: 
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = -1.2099, df = 8, p-value = 0.2608
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8198927  0.3141142
## sample estimates:
##        cor 
## -0.3932961 
## 
## Adjusted pvalue = 1

Fungi

for (plant_part in c("rhizosphere_soil", "leaves", "flowers")){
  writeLines(paste("Running", plant_part, "\n", sep=" "))
  # Subset phyloseq
  physeq.sub = subset_samples(physeqITS.fungi, Sample_type == plant_part & Field != "Freeville")
  
  # Remove fields with fewer than 3 replicates
  rep.n = data.frame(sample_data(physeq.sub), stringsAsFactors = F) %>%
    group_by(Field) %>%
    mutate(N = n()) %>%
    ungroup
  physeq.sub = prune_samples(rep.n[rep.n$N >= 3,]$sample_ID, physeq.sub)
  writeLines(paste("Removed fields:", unique(rep.n[rep.n$N < 3,]$Field)))
  
  # rarefy and normalize counts
  set.seed(72)
  physeq.sub = rarefy_even_depth(physeq.sub)
  physeq.sub = transform_sample_counts(physeq.sub, function(x) x / sum(x))
  
  # get sample metadata and split into 2 similar dataframes for contrasts
  metadata = data.frame(sample_data(physeq.sub), stringsAsFactors = F)
  meta1 = metadata %>%
    dplyr::select(sample_ID, Field)
  colnames(meta1) = c("sample1", "field1")
  meta2 = metadata %>%
    dplyr::select(sample_ID, Field)
  colnames(meta2) = c("sample2", "field2")
  
  # Calculate pairwise Bray-Curtis distances and add metadata
  BC.dist = data.frame(as.matrix(phyloseq::distance(physeq.sub, method="bray")))
  colnames(BC.dist) = gsub("X", "", gsub("\\.", "-", colnames(BC.dist)))
  dist.df = BC.dist %>%
    rownames_to_column(var="sample1") %>%
    gather(key="sample2", value="BC_dist", -sample1) %>%
    filter(sample1 != sample2) %>%
    mutate(contrast = paste(pmin(sample1, sample2), pmax(sample1, sample2), sep="+")) %>%
    separate(contrast, into=c("sample1", "sample2"), sep="\\+", remove=F) %>%
    unique() %>%
    inner_join(meta1) %>%
    inner_join(meta2) %>%
    filter(field1 != field2) %>%
    group_by(field1, field2) %>%
    summarize(mean_BC = mean(BC_dist),
              sd_BC = sd(BC_dist),
              n.con = n())
  
  # Add field data
  field.diff1 = field.data %>%
    dplyr::select(Field, Clean_Grain_Yield)
  colnames(field.diff1) = c("field1", "CGY1")
  field.diff2 = field.data %>%
    dplyr::select(Field, Clean_Grain_Yield)
  colnames(field.diff2) = c("field2", "CGY2")
  
  dist.field.df = inner_join(dist.df, field.diff1, by="field1")
  dist.field.df = inner_join(dist.field.df, field.diff2, by="field2")
  dist.field.df = dist.field.df %>%
    mutate(CGY_diff = abs(CGY1-CGY2))
  
  # Calculate correlation coeficient
  writeLines("Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield")
  corr = cor.test(dist.field.df$mean_BC, dist.field.df$CGY_diff, method = "pearson", conf.level = 0.95)
  print(corr)
  writeLines(paste("Adjusted pvalue =", p.adjust(corr$p.value, method="bonferroni", n=3), "\n"))
  
  # Plot
  title = paste(plant_part, ":", " r=", round(corr$estimate, digits=5), ", df=", corr$parameter, ", p=", round(p.adjust(corr$p.value, method="bonferroni", n=3), digits=5), sep="")
  dist.field.plot = ggplot(data=dist.field.df, aes(x=mean_BC, y=CGY_diff)) +
    geom_point() +
    geom_smooth(method="lm") +
    labs(x="Mean Bray-Curtis dissimilarity", y="Difference in clean grain yield") +
    ggtitle(label = title) +
    theme_bw()
  
  plot(dist.field.plot)
  
}
## Running rhizosphere_soil 
## 
## Removed fields: Blight
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = 0.077962, df = 4, p-value = 0.9416
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7978304  0.8244507
## sample estimates:
##        cor 
## 0.03895165 
## 
## Adjusted pvalue = 1

## Running leaves 
## 
## Removed fields: 
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = 0.78254, df = 8, p-value = 0.4564
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4362097  0.7674333
## sample estimates:
##       cor 
## 0.2666524 
## 
## Adjusted pvalue = 1

## Running flowers 
## 
## Removed fields: McGowan_Late
## Correlation analysis between mean Bray-Curtis dissimilarity and difference in clean grain yield
## 
##  Pearson's product-moment correlation
## 
## data:  dist.field.df$mean_BC and dist.field.df$CGY_diff
## t = 0.19025, df = 1, p-value = 0.8803
## alternative hypothesis: true correlation is not equal to 0
## sample estimates:
##       cor 
## 0.1868999 
## 
## Adjusted pvalue = 1

Session info

sessionInfo()
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## Random number generation:
##  RNG:     Mersenne-Twister 
##  Normal:  Inversion 
##  Sample:  Rounding 
##  
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggplot2_3.3.0               agricolae_1.3-2            
##  [3] knitr_1.28                  ape_5.3                    
##  [5] DESeq2_1.26.0               SummarizedExperiment_1.16.1
##  [7] DelayedArray_0.12.2         BiocParallel_1.20.1        
##  [9] matrixStats_0.55.0          Biobase_2.46.0             
## [11] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
## [13] IRanges_2.20.2              S4Vectors_0.24.3           
## [15] BiocGenerics_0.32.0         vegan_2.5-6                
## [17] lattice_0.20-38             permute_0.9-5              
## [19] multcomp_1.4-13             TH.data_1.0-10             
## [21] MASS_7.3-51.4               survival_3.1-8             
## [23] mvtnorm_1.1-0               kableExtra_1.1.0           
## [25] phyloseq_1.30.0             tidyr_1.0.2                
## [27] tibble_2.1.3                dplyr_0.8.4                
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.5        Hmisc_4.3-1            plyr_1.8.5            
##   [4] igraph_1.2.4.2         splines_3.6.2          AlgDesign_1.2.0       
##   [7] digest_0.6.24          foreach_1.4.8          htmltools_0.4.0       
##  [10] magrittr_1.5           checkmate_2.0.0        memoise_1.1.0         
##  [13] cluster_2.1.0          Biostrings_2.54.0      readr_1.3.1           
##  [16] annotate_1.64.0        sandwich_2.5-1         jpeg_0.1-8.1          
##  [19] colorspace_1.4-1       blob_1.2.1             rvest_0.3.5           
##  [22] xfun_0.12              libcoin_1.0-5          crayon_1.3.4          
##  [25] RCurl_1.98-1.1         jsonlite_1.6.1         genefilter_1.68.0     
##  [28] zoo_1.8-8              iterators_1.0.12       glue_1.3.1            
##  [31] gtable_0.3.0           zlibbioc_1.32.0        XVector_0.26.0        
##  [34] webshot_0.5.2          questionr_0.7.0        Rhdf5lib_1.8.0        
##  [37] dunn.test_1.3.5        scales_1.1.0           DBI_1.1.0             
##  [40] miniUI_0.1.1.1         Rcpp_1.0.3             viridisLite_0.3.0     
##  [43] xtable_1.8-4           htmlTable_1.13.3       foreign_0.8-72        
##  [46] bit_1.1-15.2           Formula_1.2-3          rcompanion_2.3.25     
##  [49] htmlwidgets_1.5.1      httr_1.4.1             RColorBrewer_1.1-2    
##  [52] modeltools_0.2-23      acepack_1.4.1          ellipsis_0.3.0        
##  [55] pkgconfig_2.0.3        XML_3.99-0.3           farver_2.0.3          
##  [58] multcompView_0.1-8     nnet_7.3-12            locfit_1.5-9.1        
##  [61] tidyselect_1.0.0       labeling_0.3           rlang_0.4.4           
##  [64] reshape2_1.4.3         later_1.0.0            AnnotationDbi_1.48.0  
##  [67] munsell_0.5.0          tools_3.6.2            RSQLite_2.2.0         
##  [70] ade4_1.7-15            EMT_1.1                evaluate_0.14         
##  [73] biomformat_1.14.0      stringr_1.4.0          fastmap_1.0.1         
##  [76] yaml_2.2.1             bit64_0.9-7            purrr_0.3.3           
##  [79] coin_1.3-1             nlme_3.1-142           mime_0.9              
##  [82] xml2_1.2.2             compiler_3.6.2         rstudioapi_0.11       
##  [85] png_0.1-7              klaR_0.6-15            geneplotter_1.64.0    
##  [88] DescTools_0.99.36      stringi_1.4.5          highr_0.8             
##  [91] Matrix_1.2-18          multtest_2.42.0        vctrs_0.2.2           
##  [94] pillar_1.4.3           lifecycle_0.1.0        lmtest_0.9-37         
##  [97] combinat_0.0-8         data.table_1.12.8      cowplot_1.0.0         
## [100] bitops_1.0-6           httpuv_1.5.2           R6_2.4.1              
## [103] latticeExtra_0.6-29    promises_1.1.0         gridExtra_2.3         
## [106] codetools_0.2-16       boot_1.3-23            assertthat_0.2.1      
## [109] rhdf5_2.30.1           nortest_1.0-4          withr_2.1.2           
## [112] GenomeInfoDbData_1.2.2 expm_0.999-4           mgcv_1.8-31           
## [115] hms_0.5.3              grid_3.6.2             rpart_4.1-15          
## [118] rmarkdown_2.1          FSA_0.8.30             shiny_1.4.0.2         
## [121] base64enc_0.1-3